/[MITgcm]/manual/s_autodiff/text/doc_ad_2.tex
ViewVC logotype

Contents of /manual/s_autodiff/text/doc_ad_2.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download) (as text)
Fri Oct 5 22:22:20 2001 UTC (23 years, 9 months ago) by heimbach
Branch: MAIN
Changes since 1.3: +292 -152 lines
File MIME type: application/x-tex
Revised part 5.

1 % $Header: /u/gcmpack/mitgcmdoc/part5/doc_ad_2.tex,v 1.3 2001/09/27 02:00:24 cnh Exp $
2 % $Name: $
3
4 {\sf Automatic differentiation} (AD), also referred to as algorithmic
5 (or, more loosely, computational) differentiation, involves
6 automatically deriving code to calculate
7 partial derivatives from an existing fully non-linear prognostic code.
8 (see \cite{gri:00}).
9 A software tool is used that parses and transforms source files
10 according to a set of linguistic and mathematical rules.
11 AD tools are like source-to-source translators in that
12 they parse a program code as input and produce a new program code
13 as output.
14 However, unlike a pure source-to-source translation, the output program
15 represents a new algorithm, such as the evaluation of the
16 Jacobian, the Hessian, or higher derivative operators.
17 In principle, a variety of derived algorithms
18 can be generated automatically in this way.
19
20 The MITGCM has been adapted for use with the
21 Tangent linear and Adjoint Model Compiler (TAMC) and its successor TAF
22 (Transformation of Algorithms in Fortran), developed
23 by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}).
24 The first application of the adjoint of the MITGCM for senistivity
25 studies has been published by \cite{maro-eta:99}.
26 \cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint
27 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
32 derivative of a function with
33 respect to a set of input variables.
34 Treating a given forward code as a composition of operations --
35 each line representing a compositional element, the chain rule is
36 rigorously applied to the code, line by line. The resulting
37 tangent linear or adjoint code,
38 then, may be thought of as the composition in
39 forward or reverse order, respectively, of the
40 Jacobian matrices of the forward code's compositional elements.
41
42 %**********************************************************************
43 \section{Some basic algebra}
44 \label{sec_ad_algebra}
45 %**********************************************************************
46
47 Let $ \cal{M} $ be a general nonlinear, model, i.e. a
48 mapping from the $m$-dimensional space
49 $U \subset I\!\!R^m$ of input variables
50 $\vec{u}=(u_1,\ldots,u_m)$
51 (model parameters, initial conditions, boundary conditions
52 such as forcing functions) to the $n$-dimensional space
53 $V \subset I\!\!R^n$ of
54 model output variable $\vec{v}=(v_1,\ldots,v_n)$
55 (model state, model diagnostcs, objective function, ...)
56 under consideration,
57 %
58 \begin{equation}
59 \begin{split}
60 {\cal M} \, : & \, U \,\, \longrightarrow \, V \\
61 ~ & \, \vec{u} \,\, \longmapsto \, \vec{v} \, = \,
62 {\cal M}(\vec{u})
63 \label{fulloperator}
64 \end{split}
65 \end{equation}
66 %
67 The vectors $ \vec{u} \in U $ and $ v \in V $ may be represented w.r.t.
68 some given basis vectors
69 $ {\rm span} (U) = \{ {\vec{e}_i} \}_{i = 1, \ldots , m} $ and
70 $ {\rm span} (V) = \{ {\vec{f}_j} \}_{j = 1, \ldots , n} $ as
71 \[
72 \vec{u} \, = \, \sum_{i=1}^{m} u_i \, {\vec{e}_i},
73 \qquad
74 \vec{v} \, = \, \sum_{j=1}^{n} v_j \, {\vec{f}_j}
75 \]
76
77 Two routes may be followed to determine the sensitivity of the
78 output variable $\vec{v}$ to its input $\vec{u}$.
79
80 \subsection{Forward or direct sensitivity}
81 %
82 Consider a perturbation to the input variables $\delta \vec{u}$
83 (typically a single component
84 $\delta \vec{u} = \delta u_{i} \, {\vec{e}_{i}}$).
85 Their effect on the output may be obtained via the linear
86 approximation of the model $ {\cal M}$ in terms of its Jacobian matrix
87 $ M $, evaluated in the point $u^{(0)}$ according to
88 %
89 \begin{equation}
90 \delta \vec{v} \, = \, M |_{\vec{u}^{(0)}} \, \delta \vec{u}
91 \label{tangent_linear}
92 \end{equation}
93 with resulting output perturbation $\delta \vec{v}$.
94 In components
95 $M_{j i} \, = \, \partial {\cal M}_{j} / \partial u_{i} $,
96 it reads
97 %
98 \begin{equation}
99 \delta v_{j} \, = \, \sum_{i}
100 \left. \frac{\partial {\cal M}_{j}}{\partial u_{i}} \right|_{u^{(0)}} \,
101 \delta u_{i}
102 \label{jacobi_matrix}
103 \end{equation}
104 %
105 Eq. (\ref{tangent_linear}) is the {\sf tangent linear model (TLM)}.
106 In contrast to the full nonlinear model $ {\cal M} $, the operator
107 $ M $ is just a matrix
108 which can readily be used to find the forward sensitivity of $\vec{v}$ to
109 perturbations in $u$,
110 but if there are very many input variables $(\gg O(10^{6})$ for
111 large-scale oceanographic application), it quickly becomes
112 prohibitive to proceed directly as in (\ref{tangent_linear}),
113 if the impact of each component $ {\bf e_{i}} $ is to be assessed.
114
115 \subsection{Reverse or adjoint sensitivity}
116 %
117 Let us consider the special case of a
118 scalar objective function ${\cal J}(\vec{v})$ of the model output (e.g.
119 the total meridional heat transport,
120 the total uptake of $CO_{2}$ in the Southern
121 Ocean over a time interval,
122 or a measure of some model-to-data misfit)
123 %
124 \begin{eqnarray}
125 \begin{array}{cccccc}
126 {\cal J} \, : & U &
127 \longrightarrow & V &
128 \longrightarrow & I \!\! R \\
129 ~ & \vec{u} & \longmapsto & \vec{v}={\cal M}(\vec{u}) &
130 \longmapsto & {\cal J}(\vec{u}) = {\cal J}({\cal M}(\vec{u}))
131 \end{array}
132 \label{compo}
133 \end{eqnarray}
134 %
135 The perturbation of $ {\cal J} $ around a fixed point $ {\cal J}_0 $,
136 \[
137 {\cal J} \, = \, {\cal J}_0 \, + \, \delta {\cal J}
138 \]
139 can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $
140 w.r.t. their corresponding inner product
141 $\left\langle \,\, , \,\, \right\rangle $
142 %
143 \begin{equation}
144 \begin{split}
145 {\cal J} & = \,
146 {\cal J} |_{\vec{u}^{(0)}} \, + \,
147 \left\langle \, \nabla _{u}{\cal J}^T |_{\vec{u}^{(0)}} \, , \, \delta \vec{u} \, \right\rangle
148 \, + \, O(\delta \vec{u}^2) \\
149 ~ & = \,
150 {\cal J} |_{\vec{v}^{(0)}} \, + \,
151 \left\langle \, \nabla _{v}{\cal J}^T |_{\vec{v}^{(0)}} \, , \, \delta \vec{v} \, \right\rangle
152 \, + \, O(\delta \vec{v}^2)
153 \end{split}
154 \label{deljidentity}
155 \end{equation}
156 %
157 (note, that the gradient $ \nabla f $ is a co-vector, therefore
158 its transpose is required in the above inner product).
159 Then, using the representation of
160 $ \delta {\cal J} =
161 \left\langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \right\rangle $,
162 the definition
163 of an adjoint operator $ A^{\ast} $ of a given operator $ A $,
164 \[
165 \left\langle \, A^{\ast} \vec{x} \, , \, \vec{y} \, \right\rangle =
166 \left\langle \, \vec{x} \, , \, A \vec{y} \, \right\rangle
167 \]
168 which for finite-dimensional vector spaces is just the
169 transpose of $ A $,
170 \[
171 A^{\ast} \, = \, A^T
172 \]
173 and from eq. (\ref{tangent_linear}), (\ref{deljidentity}),
174 we note that
175 (omitting $|$'s):
176 %
177 \begin{equation}
178 \delta {\cal J}
179 \, = \,
180 \left\langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \right\rangle
181 \, = \,
182 \left\langle \, \nabla _{v}{\cal J}^T \, , \, M \, \delta \vec{u} \, \right\rangle
183 \, = \,
184 \left\langle \, M^T \, \nabla _{v}{\cal J}^T \, , \,
185 \delta \vec{u} \, \right\rangle
186 \label{inner}
187 \end{equation}
188 %
189 With the identity (\ref{deljidentity}), we then find that
190 the gradient $ \nabla _{u}{\cal J} $ can be readily inferred by
191 invoking the adjoint $ M^{\ast } $ of the tangent linear model $ M $
192 %
193 \begin{equation}
194 \begin{split}
195 \nabla _{u}{\cal J}^T |_{\vec{u}} &
196 = \, M^T |_{\vec{u}} \cdot \nabla _{v}{\cal J}^T |_{\vec{v}} \\
197 ~ & = \, M^T |_{\vec{u}} \cdot \delta \vec{v}^{\ast} \\
198 ~ & = \, \delta \vec{u}^{\ast}
199 \end{split}
200 \label{adjoint}
201 \end{equation}
202 %
203 Eq. (\ref{adjoint}) is the {\sf adjoint model (ADM)},
204 in which $M^T$ is the adjoint (here, the transpose) of the
205 tangent linear operator $M$, $ \delta \vec{v}^{\ast} $
206 the adjoint variable of the model state $ \vec{v} $, and
207 $ \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
210 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 funciton 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}
230 {\cal J}({\cal M}(\vec{u})) \, = \,
231 {\cal J} ( {\cal M}_{\Lambda} ( {\cal M}_{\Lambda-1} (
232 ...... ( {\cal M}_{\lambda} (
233 ......
234 ( {\cal M}_{1} ( {\cal M}_{0}(\vec{u}) )))))
235 \label{compos}
236 \end{equation}
237 %
238 Then, according to the chain rule, the forward calculation reads,
239 in terms of the Jacobi matrices
240 (we've omitted the $ | $'s which, nevertheless are important
241 to the aspect of {\it tangent} linearity;
242 note also that by definition
243 $ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle
244 = \nabla_v {\cal J} \cdot \delta \vec{v} $ )
245 %
246 \begin{equation}
247 \begin{split}
248 \nabla_v {\cal J} (M(\delta \vec{u})) & = \,
249 \nabla_v {\cal J} \cdot M_{\Lambda}
250 \cdot ...... \cdot M_{\lambda} \cdot ...... \cdot
251 M_{1} \cdot M_{0} \cdot \delta \vec{u} \\
252 ~ & = \, \nabla_v {\cal J} \cdot \delta \vec{v} \\
253 \end{split}
254 \label{forward}
255 \end{equation}
256 %
257 whereas in reverse mode we have
258 %
259 \begin{equation}
260 \boxed{
261 \begin{split}
262 M^T ( \nabla_v {\cal J}^T) & = \,
263 M_{0}^T \cdot M_{1}^T
264 \cdot ...... \cdot M_{\lambda}^T \cdot ...... \cdot
265 M_{\Lambda}^T \cdot \nabla_v {\cal J}^T \\
266 ~ & = \, M_{0}^T \cdot M_{1}^T
267 \cdot ...... \cdot
268 \nabla_{v^{(\lambda)}} {\cal J}^T \\
269 ~ & = \, \nabla_u {\cal J}^T
270 \end{split}
271 }
272 \label{reverse}
273 \end{equation}
274 %
275 clearly expressing the reverse nature of the calculation.
276 Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers.
277 If the intermediate steps $\lambda$ in
278 eqn. (\ref{compos}) -- (\ref{reverse})
279 represent the model state (forward or adjoint) at each
280 intermediate time step as noted above, then correspondingly,
281 $ M^T (\delta \vec{v}^{(\lambda) \, \ast}) =
282 \delta \vec{v}^{(\lambda-1) \, \ast} $ for the adjoint variables.
283 It thus becomes evident that the adjoint calculation also
284 yields the adjoint of each model state component
285 $ \vec{v}^{(\lambda)} $ at each intermediate step $ \lambda $, namely
286 %
287 \begin{equation}
288 \boxed{
289 \begin{split}
290 \nabla_{v^{(\lambda)}} {\cal J}^T |_{\vec{v}^{(\lambda)}}
291 & = \,
292 M_{\lambda}^T |_{\vec{v}^{(\lambda)}} \cdot ...... \cdot
293 M_{\Lambda}^T |_{\vec{v}^{(\lambda)}} \cdot \delta \vec{v}^{\ast} \\
294 ~ & = \, \delta \vec{v}^{(\lambda) \, \ast}
295 \end{split}
296 }
297 \end{equation}
298 %
299 in close analogy to eq. (\ref{adjoint})
300 We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$
301 are the Lagrange multipliers of the model equations which determine
302 $ \vec{v}^{(\lambda)}$.
303
304 In coponents, eq. (\ref{adjoint}) reads as follows.
305 Let
306 \[
307 \begin{array}{rclcrcl}
308 \delta \vec{u} & = &
309 \left( \delta u_1,\ldots, \delta u_m \right)^T , & \qquad &
310 \delta \vec{u}^{\ast} \,\, = \,\, \nabla_u {\cal J}^T & = &
311 \left(
312 \frac{\partial {\cal J}}{\partial u_1},\ldots,
313 \frac{\partial {\cal J}}{\partial u_m}
314 \right)^T \\
315 \delta \vec{v} & = &
316 \left( \delta v_1,\ldots, \delta u_n \right)^T , & \qquad &
317 \delta \vec{v}^{\ast} \,\, = \,\, \nabla_v {\cal J}^T & = &
318 \left(
319 \frac{\partial {\cal J}}{\partial v_1},\ldots,
320 \frac{\partial {\cal J}}{\partial v_n}
321 \right)^T \\
322 \end{array}
323 \]
324 denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,
325 and their adjoint varaiables;
326 further
327 \[
328 M \, = \, \left(
329 \begin{array}{ccc}
330 \frac{\partial {\cal M}_1}{\partial u_1} & \ldots &
331 \frac{\partial {\cal M}_1}{\partial u_m} \\
332 \vdots & ~ & \vdots \\
333 \frac{\partial {\cal M}_n}{\partial u_1} & \ldots &
334 \frac{\partial {\cal M}_n}{\partial u_m} \\
335 \end{array}
336 \right)
337 \]
338 is the Jacobi matrix of $ {\cal M} $
339 (an $ n \times m $ matrix)
340 such that $ \delta \vec{v} = M \cdot \delta \vec{u} $, or
341 \[
342 \delta v_{j}
343 \, = \, \sum_{i=1}^m M_{ji} \, \delta u_{i}
344 \, = \, \sum_{i=1}^m \, \frac{\partial {\cal M}_{j}}{\partial u_{i}}
345 \delta u_{i}
346 \]
347 %
348 Then eq. (\ref{adjoint}) takes the form
349 \[
350 \delta u_{i}^{\ast}
351 \, = \, \sum_{j=1}^n M_{ji} \, \delta v_{j}^{\ast}
352 \, = \, \sum_{j=1}^n \, \frac{\partial {\cal M}_{j}}{\partial u_{i}}
353 \delta v_{j}^{\ast}
354 \]
355 %
356 or
357 %
358 \[
359 \left(
360 \begin{array}{c}
361 \left. \frac{\partial}{\partial u_1} {\cal J} \right|_{\vec{u}^{(0)}} \\
362 \vdots \\
363 \left. \frac{\partial}{\partial u_m} {\cal J} \right|_{\vec{u}^{(0)}} \\
364 \end{array}
365 \right)
366 \, = \,
367 \left(
368 \begin{array}{ccc}
369 \left. \frac{\partial {\cal M}_1}{\partial u_1} \right|_{\vec{u}^{(0)}}
370 & \ldots &
371 \left. \frac{\partial {\cal M}_n}{\partial u_1} \right|_{\vec{u}^{(0)}} \\
372 \vdots & ~ & \vdots \\
373 \left. \frac{\partial {\cal M}_1}{\partial u_m} \right|_{\vec{u}^{(0)}}
374 & \ldots &
375 \left. \frac{\partial {\cal M}_n}{\partial u_m} \right|_{\vec{u}^{(0)}} \\
376 \end{array}
377 \right)
378 \cdot
379 \left(
380 \begin{array}{c}
381 \left. \frac{\partial}{\partial v_1} {\cal J} \right|_{\vec{v}} \\
382 \vdots \\
383 \left. \frac{\partial}{\partial v_n} {\cal J} \right|_{\vec{v}} \\
384 \end{array}
385 \right)
386 \]
387 %
388 Furthermore, the adjoint $ \delta v^{(\lambda) \, \ast} $
389 of any intermediate state $ v^{(\lambda)} $
390 may be obtained, using the intermediate Jacobian
391 (an $ n_{\lambda+1} \times n_{\lambda} $ matrix)
392 %
393 \[
394 M_{\lambda} \, = \,
395 \left(
396 \begin{array}{ccc}
397 \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}
398 & \ldots &
399 \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}} \\
400 \vdots & ~ & \vdots \\
401 \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1}
402 & \ldots &
403 \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\
404 \end{array}
405 \right)
406 \]
407 %
408 and the shorthand notation for the adjoint variables
409 $ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}}
410 {\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $,
411 for intermediate components, yielding
412 \begin{equation}
413 \small
414 \begin{split}
415 \left(
416 \begin{array}{c}
417 \delta v^{(\lambda) \, \ast}_1 \\
418 \vdots \\
419 \delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\
420 \end{array}
421 \right)
422 \, = &
423 \left(
424 \begin{array}{ccc}
425 \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}
426 & \ldots \,\, \ldots &
427 \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\
428 \vdots & ~ & \vdots \\
429 \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}}
430 & \ldots \,\, \ldots &
431 \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\
432 \end{array}
433 \right)
434 \cdot
435 %
436 \\ ~ & ~
437 \\ ~ &
438 %
439 \left(
440 \begin{array}{ccc}
441 \frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_1}
442 & \ldots &
443 \frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_1} \\
444 \vdots & ~ & \vdots \\
445 \vdots & ~ & \vdots \\
446 \frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}}
447 & \ldots &
448 \frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}} \\
449 \end{array}
450 \right)
451 \cdot \, \ldots \, \cdot
452 \left(
453 \begin{array}{c}
454 \delta v^{\ast}_1 \\
455 \vdots \\
456 \delta v^{\ast}_{n} \\
457 \end{array}
458 \right)
459 \end{split}
460 \end{equation}
461
462 Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in
463 showing the advantage of the reverse over the forward mode
464 if the gradient $\nabla _{u}{\cal J}$, i.e. the sensitivity of the
465 cost function $ {\cal J} $ with respect to {\it all} input
466 variables $u$
467 (or the sensitivity of the cost function with respect to
468 {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.
469 In order to be able to solve for each component of the gradient
470 $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})
471 a forward calulation has to be performed for each component seperately,
472 i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $
473 for the $i$-th forward calculation.
474 Then, (\ref{forward}) represents the
475 projection of $ \nabla_u {\cal J} $ onto the $i$-th component.
476 The full gradient is retrieved from the $ m $ forward calculations.
477 In contrast, eq. (\ref{reverse}) yields the full
478 gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients
479 $\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation.
480
481 Note, that if $ {\cal J} $ is a vector-valued function
482 of dimension $ l > 1 $,
483 eq. (\ref{reverse}) has to be modified according to
484 \[
485 M^T \left( \nabla_v {\cal J}^T \left(\delta \vec{J}\right) \right)
486 \, = \,
487 \nabla_u {\cal J}^T \cdot \delta \vec{J}
488 \]
489 where now $ \delta \vec{J} \in I\!\!R^l $ is a vector of
490 dimenison $ l $.
491 In this case $ l $ reverse simulations have to be performed
492 for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.
493 Then, the reverse mode is more efficient as long as
494 $ l < n $, otherwise the forward mode is preferable.
495 Stricly, the reverse mode is called adjoint mode only for
496 $ l = 1 $.
497
498 A detailed analysis of the underlying numerical operations
499 shows that the computation of $\nabla _{u}{\cal J}$ in this way
500 requires about 2 to 5 times the computation of the cost function.
501 Alternatively, the gradient vector could be approximated
502 by finite differences, requiring $m$ computations
503 of the perturbed cost function.
504
505 To conclude we give two examples of commonly used types
506 of cost functions:
507
508 \paragraph{Example 1:
509 $ {\cal J} = v_{j} (T) $} ~ \\
510 The cost function consists of the $j$-th component of the model state
511 $ \vec{v} $ at time $T$.
512 Then $ \nabla_v {\cal J}^T = {\vec{f}_{j}} $ is just the $j$-th
513 unit vector. The $ \nabla_u {\cal J}^T $
514 is the projection of the adjoint
515 operator onto the $j$-th component ${\bf f_{j}}$,
516 \[
517 \nabla_u {\cal J}^T
518 \, = \, M^T \cdot \nabla_v {\cal J}^T
519 \, = \, \sum_{i} M^T_{ji} \, {\vec{e}_{i}}
520 \]
521
522 \paragraph{Example 2:
523 $ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, ,
524 \, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\
525 The cost function represents the quadratic model vs. data misfit.
526 Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the
527 operator which maps the model state space onto the data space.
528 Then, $ \nabla_v {\cal J} $ takes the form
529 %
530 \begin{equation*}
531 \begin{split}
532 \nabla_v {\cal J}^T & = \, 2 \, \, H \cdot
533 \left( \, {\cal H}(\vec{v}) - \vec{d} \, \right) \\
534 ~ & = \, 2 \sum_{j} \left\{ \sum_k
535 \frac{\partial {\cal H}_k}{\partial v_{j}}
536 \left( {\cal H}_k (\vec{v}) - d_k \right)
537 \right\} \, {\vec{f}_{j}} \\
538 \end{split}
539 \end{equation*}
540 %
541 where $H_{kj} = \partial {\cal H}_k / \partial v_{j} $ is the
542 Jacobi matrix of the data projection operator.
543 Thus, the gradient $ \nabla_u {\cal J} $ is given by the
544 adjoint operator,
545 driven by the model vs. data misfit:
546 \[
547 \nabla_u {\cal J}^T \, = \, 2 \, M^T \cdot
548 H \cdot \left( {\cal H}(\vec{v}) - \vec{d} \, \right)
549 \]
550
551 \subsection{Storing vs. recomputation in reverse mode}
552 \label{checkpointing}
553
554 We note an important aspect of the forward vs. reverse
555 mode calculation.
556 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
559 $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$
560 are needed to evaluate the intermediate Jacobian
561 $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.
562 In the forward mode, the intermediate results are required
563 in the same order as computed by the full forward model ${\cal M}$,
564 but in the reverse mode they are required in the reverse order.
565 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
567 calculation. Alternatively, the complete model state up to the
568 point of evaluation has to be recomputed whenever its value is required.
569
570 A method to balance the amount of recomputations vs.
571 storage requirements is called {\sf checkpointing}
572 (e.g. \cite{res-eta:98}).
573 It is depicted in \reffig{3levelcheck} for a 3-level checkpointing
574 [as an example, we give explicit numbers for a 3-day
575 integration with a 1-hourly timestep in square brackets].
576 \begin{itemize}
577 %
578 \item [$lev3$]
579 In a first step, the model trajectory is subdivided into
580 $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],
581 with the label $lev3$ for this outermost loop.
582 The model is then integrated along the full trajectory,
583 and the model state stored only at every $ k_{i}^{lev3} $-th timestep
584 [i.e. 3 times, at
585 $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].
586 %
587 \item [$lev2$]
588 In a second step each subsection itself is divided into
589 $ {n}^{lev2} $ sub-subsections
590 [$ {n}^{lev2} $=4 6-hour intervals per subsection].
591 The model picks up at the last outermost dumped state
592 $ v_{k_{n}^{lev3}} $ and is integrated forward in time along
593 the last subsection, with the label $lev2$ for this
594 intermediate loop.
595 The model state is now stored at every $ k_{i}^{lev2} $-th
596 timestep
597 [i.e. 4 times, at
598 $ i = 0,1,2,3 $ corresponding to $ k_{i}^{lev2} = 48, 54, 60, 66 $].
599 %
600 \item [$lev1$]
601 Finally, the model picks up at the last intermediate dump state
602 $ v_{k_{n}^{lev2}} $ and is integrated forward in time along
603 the last sub-subsection, with the label $lev1$ for this
604 intermediate loop.
605 Within this sub-subsection only, the model state is stored
606 at every timestep
607 [i.e. every hour $ i=0,...,5$ corresponding to
608 $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].
609 Thus, the final state $ v_n = v_{k_{n}^{lev1}} $ is reached
610 and the model state of all peceeding timesteps along the last
611 sub-subsections are available, enabling integration backwards
612 in time along the last sub-subsection.
613 Thus, the adjoint can be computed along this last
614 sub-subsection $k_{n}^{lev2}$.
615 %
616 \end{itemize}
617 %
618 This procedure is repeated consecutively for each previous
619 sub-subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $
620 carrying the adjoint computation to the initial time
621 of the subsection $k_{n}^{lev3}$.
622 Then, the procedure is repeated for the previous subsection
623 $k_{n-1}^{lev3}$
624 carrying the adjoint computation to the initial time
625 $k_{1}^{lev3}$.
626
627 For the full model trajectory of
628 $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps
629 the required storing of the model state was significantly reduced to
630 $ n^{lev1} + n^{lev2} + n^{lev3} $
631 [i.e. for the 3-day integration with a total oof 72 timesteps
632 the model state was stored 13 times].
633 This saving in memory comes at a cost of a required
634 3 full forward integrations of the model (one for each
635 checkpointing level).
636 The balance of storage vs. recomputation certainly depends
637 on the computing resources available.
638
639 \begin{figure}[t!]
640 \centering
641 %\psdraft
642 \psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}
643 \psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}
644 \psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}
645 \psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}
646 \psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}
647 \psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}
648 \psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}
649 \psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}
650 \mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}}
651 %\psfull
652 \caption
653 {Schematic view of intermediate dump and restart for
654 3-level checkpointing.}
655 \label{fig:3levelcheck}
656 \end{figure}
657
658 % \subsection{Optimal perturbations}
659 % \label{sec_optpert}
660
661
662 % \subsection{Error covariance estimate and Hessian matrix}
663 % \label{sec_hessian}
664
665 \newpage
666
667 %**********************************************************************
668 \section{AD-specific setup by example: sensitivity of carbon sequestration}
669 \label{sec_ad_setup_ex}
670 %**********************************************************************
671
672 The MITGCM has been adapted to enable AD using TAMC or TAF.
673 The present description, therefore, is specific to the
674 use of TAMC or TAF as AD tool.
675 The following sections describe the steps which are necessary to
676 generate a tangent linear or adjoint model of the MITGCM.
677 We take as an example the sensitivity of carbon sequestration
678 in the ocean.
679 The AD-relevant hooks in the code are sketched in
680 \reffig{adthemodel}, \reffig{adthemain}.
681
682 \subsection{Overview of the experiment}
683
684 We describe an adjoint sensitivity analysis of outgassing from
685 the ocean into the atmosphere of a carbon-like tracer injected
686 into the ocean interior (see \cite{hil-eta:01}).
687
688 \subsubsection{Passive tracer equation}
689
690 For this work the MITGCM was augmented with a thermodynamically
691 inactive tracer, $C$. Tracer residing in the ocean
692 model surface layer is outgassed according to a relaxation time scale,
693 $\mu$. Within the ocean interior, the tracer is passively advected
694 by the ocean model currents. The full equation for the time evolution
695 %
696 \begin{equation}
697 \label{carbon_ddt}
698 \frac{\partial C}{\partial t} \, = \,
699 -U\cdot \nabla C \, - \, \mu C \, + \, \Gamma(C) \,+ \, S
700 \end{equation}
701 %
702 also includes a source term $S$. This term
703 represents interior sources of $C$ such as would arise due to
704 direct injection.
705 The velocity term, $U$, is the sum of the
706 model Eulerian circulation and an eddy-induced velocity, the latter
707 parameterized according to Gent/McWilliams
708 (\cite{gen-mcw:90, gen-eta:95}).
709 The convection function, $\Gamma$, mixes $C$ vertically wherever the
710 fluid is locally statically unstable.
711
712 The outgassing time scale, $\mu$, in eqn. (\ref{carbon_ddt})
713 is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface
714 ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt})
715 is valid as a prognostic equation for small perturbations in oceanic
716 carbon concentrations. This configuration provides a
717 powerful tool for examining the impact of large-scale ocean circulation
718 on $ CO_2 $ outgassing due to interior injections.
719 As source we choose a constant in time injection of
720 $ S = 1 \,\, {\rm mol / s}$.
721
722 \subsubsection{Model configuration}
723
724 The model configuration employed has a constant
725 $4^\circ \times 4^\circ$ resolution horizontal grid and realistic
726 geography and bathymetry. Twenty vertical layers are used with
727 vertical spacing ranging
728 from 50 m near the surface to 815 m at depth.
729 Driven to steady-state by climatalogical wind-stress, heat and
730 fresh-water forcing the model reproduces well known large-scale
731 features of the ocean general circulation.
732
733 \subsubsection{Outgassing cost function}
734
735 To quantify and understand outgassing due to injections of $C$
736 in eqn. (\ref{carbon_ddt}),
737 we define a cost function $ {\cal J} $ that measures the total amount of
738 tracer outgassed at each timestep:
739 %
740 \begin{equation}
741 \label{cost_tracer}
742 {\cal J}(t=T)=\int_{t=0}^{t=T}\int_{A} \mu C \, dA \, dt
743 \end{equation}
744 %
745 Equation(\ref{cost_tracer}) integrates the outgassing term, $\mu C$,
746 from (\ref{carbon_ddt})
747 over the entire ocean surface area, $A$, and accumulates it
748 up to time $T$.
749 Physically, ${\cal J}$ can be thought of as representing the amount of
750 $CO_2$ that our model predicts would be outgassed following an
751 injection at rate $S$.
752 The sensitivity of ${\cal J}$ to the spatial location of $S$,
753 $\frac{\partial {\cal J}}{\partial S}$,
754 can be used to identify regions from which circulation
755 would cause $CO_2$ to rapidly outgas following injection
756 and regions in which $CO_2$ injections would remain effectively
757 sequesterd within the ocean.
758
759 \subsection{Code configuration}
760
761 The model configuration for this experiment resides under the
762 directory {\it verification/carbon/}.
763 The code customisation routines are in {\it verification/carbon/code/}:
764 %
765 \begin{itemize}
766 %
767 \item {\it .genmakerc}
768 %
769 \item {\it COST\_CPPOPTIONS.h}
770 %
771 \item {\it CPP\_EEOPTIONS.h}
772 %
773 \item {\it CPP\_OPTIONS.h}
774 %
775 \item {\it CTRL\_OPTIONS.h}
776 %
777 \item {\it ECCO\_OPTIONS.h}
778 %
779 \item {\it SIZE.h}
780 %
781 \item {\it adcommon.h}
782 %
783 \item {\it tamc.h}
784 %
785 \end{itemize}
786 %
787 The runtime flag and parameters settings are contained in
788 {\it verification/carbon/input/},
789 together with the forcing fields and and restart files:
790 %
791 \begin{itemize}
792 %
793 \item {\it data}
794 %
795 \item {\it data.cost}
796 %
797 \item {\it data.ctrl}
798 %
799 \item {\it data.gmredi}
800 %
801 \item {\it data.grdchk}
802 %
803 \item {\it data.optim}
804 %
805 \item {\it data.pkg}
806 %
807 \item {\it eedata}
808 %
809 \item {\it topog.bin}
810 %
811 \item {\it windx.bin, windy.bin}
812 %
813 \item {\it salt.bin, theta.bin}
814 %
815 \item {\it SSS.bin, SST.bin}
816 %
817 \item {\it pickup*}
818 %
819 \end{itemize}
820 %
821 Finally, the file to generate the adjoint code resides in
822 $ adjoint/ $:
823 %
824 \begin{itemize}
825 %
826 \item {\it makefile}
827 %
828 \end{itemize}
829 %
830
831 Below we describe the customisations of this files which are
832 specific to this experiment.
833
834 \subsubsection{File {\it .genmakerc}}
835 This file overwrites default settings of {\it genmake}.
836 In the present example it is used to switch on the following
837 packages which are related to automatic differentiation
838 and are disabled by default: \\
839 \hspace*{4ex} {\tt set ENABLE=( autodiff cost ctrl ecco gmredi grdchk kpp )} \\
840 Other packages which are not needed are switched off: \\
841 \hspace*{4ex} {\tt set DISABLE=( aim obcs zonal\_filt shap\_filt cal exf )}
842
843 \subsubsection{File {\it COST\_CPPOPTIONS.h, CTRL\_OPTIONS.h}}
844
845 These files used to contain package-specific CPP-options
846 (see Section \ref{???}).
847 For technical reasons those options have been grouped together
848 in the file {\it ECCO\_OPTIONS.h}.
849 To retain the modularity, the files have been kept and contain
850 the standard include of the {\it CPP\_OPTIONS.h} file.
851
852 \subsubsection{File {\it CPP\_EEOPTIONS.h}}
853
854 This file contains 'wrapper'-specific CPP options.
855 It only needs to be changed if the code is to be run
856 in a parallel environment (see Section \ref{???}).
857
858 \subsubsection{File {\it CPP\_OPTIONS.h}}
859
860 This file contains model-specific CPP options
861 (see Section \ref{???}).
862 Most options are related to the forward model setup.
863 They are identical to the global steady circulation setup of
864 {\it verification/exp2/}.
865 The three options specific to this experiment are \\
866 \hspace*{4ex} {\tt \#define ALLOW\_PASSIVE\_TRACER} \\
867 This flag enables the code to carry through the
868 advection/diffusion of a passive tracer along the
869 model integration. \\
870 \hspace*{4ex} {\tt \#define ALLOW\_MIT\_ADJOINT\_RUN} \\
871 This flag enables the inclusion of some AD-related fields
872 concerning initialisation, link between control variables
873 and forward model variables, and the call to the top-level
874 forward/adjoint subroutine {\it adthe\_main\_loop}
875 instead of {\it the\_main\_loop}. \\
876 \hspace*{4ex} {\tt \#define ALLOW\_GRADIENT\_CHECK} \\
877 This flag enables the gradient check package.
878 After computing the unperturbed cost function and its gradient,
879 a series of computations are performed for which \\
880 $\bullet$ an element of the control vector is perturbed \\
881 $\bullet$ the cost function w.r.t. the perturbed element is
882 computed \\
883 $\bullet$ the difference between the perturbed and unperturbed
884 cost function is computed to compute the finite difference gradient \\
885 $\bullet$ the finite difference gradient is compared with the
886 adjoint-generated gradient.
887 The gradient check package is further described in Section ???.
888
889 \subsubsection{File {\it ECCO\_OPTIONS.h}}
890
891 The CPP options of several AD-related packages are grouped
892 in this file:
893 %
894 \begin{itemize}
895 %
896 \item
897 Adjoint support package: {\it pkg/autodiff/} \\
898 This package contains hand-written adjoint code such as
899 active file handling, flow directives for files which must not
900 be differentiated, and TAMC-specific header files. \\
901 \hspace*{4ex} {\tt \#define ALLOW\_AUTODIFF\_TAMC} \\
902 defines TAMC-related features in the code. \\
903 \hspace*{4ex} {\tt \#define ALLOW\_TAMC\_CHECKPOINTING} \\
904 enables the checkpointing feature of TAMC
905 (see Section \ref{???}).
906 In the present example a 3-level checkpointing is implemented.
907 The code contains the relevant store directives, common block
908 and tape initialisations, storing key computation,
909 and loop index handling.
910 The checkpointing length at each level is defined in
911 file {\it tamc.h}, cf. below.
912 %
913 \item Cost function package: {\it pkg/cost/} \\
914 This package contains all relevant routines for
915 initialising, accumulating and finalizing the cost function
916 (see Section \ref{???}). \\
917 \hspace*{4ex} {\tt \#define ALLOW\_COST} \\
918 enables all general aspects of the cost function handling,
919 in particular the hooks in the foorward code for
920 initialising, accumulating and finalizing the cost function. \\
921 \hspace*{4ex} {\tt \#define ALLOW\_COST\_TRACER} \\
922 includes the call to the cost function for this
923 particular experiment, eqn. (\ref{cost_tracer}).
924 %
925 \item Control variable package: {\it pkg/ctrl/} \\
926 This package contains all relevant routines for
927 the handling of the control vector.
928 Each control variable can be enabled/disabled with its own flag: \\
929 \begin{tabular}{ll}
930 \hspace*{2ex} {\tt \#define ALLOW\_THETA0\_CONTROL} &
931 initial temperature \\
932 \hspace*{2ex} {\tt \#define ALLOW\_SALT0\_CONTROL} &
933 initial salinity \\
934 \hspace*{2ex} {\tt \#define ALLOW\_TR0\_CONTROL} &
935 initial passive tracer concentration \\
936 \hspace*{2ex} {\tt \#define ALLOW\_TAUU0\_CONTROL} &
937 zonal wind stress \\
938 \hspace*{2ex} {\tt \#define ALLOW\_TAUV0\_CONTROL} &
939 meridional wind stress \\
940 \hspace*{2ex} {\tt \#define ALLOW\_SFLUX0\_CONTROL} &
941 freshwater flux \\
942 \hspace*{2ex} {\tt \#define ALLOW\_HFLUX0\_CONTROL} &
943 heat flux \\
944 \hspace*{2ex} {\tt \#define ALLOW\_DIFFKR\_CONTROL} &
945 diapycnal diffusivity \\
946 \hspace*{2ex} {\tt \#undef ALLOW\_KAPPAGM\_CONTROL} &
947 isopycnal diffusivity \\
948 \end{tabular}
949 %
950 \end{itemize}
951
952 \subsubsection{File {\it SIZE.h}}
953
954 The file contains the grid point dimensions of the forward
955 model. It is identical to the {\it verification/exp2/}: \\
956 \hspace*{4ex} {\tt sNx = 90} \\
957 \hspace*{4ex} {\tt sNy = 40} \\
958 \hspace*{4ex} {\tt Nr = 20} \\
959 It correpsponds to a single-tile/single-processor setup:
960 {\tt nSx = nSy = 1, nPx = nPy = 1},
961 with standard overlap dimensioning
962 {\tt OLx = OLy = 3}.
963
964 \subsubsection{File {\it adcommon.h}}
965
966 This file contains common blocks of some adjoint variables
967 that are generated by TAMC.
968 The common blocks are used by the adjoint support routine
969 {\it addummy\_in\_stepping} which needs to access those variables:
970
971 \begin{tabular}{ll}
972 \hspace*{4ex} {\tt common /addynvars\_r/} &
973 \hspace*{4ex} is related to {\it DYNVARS.h} \\
974 \hspace*{4ex} {\tt common /addynvars\_cd/} &
975 \hspace*{4ex} is related to {\it DYNVARS.h} \\
976 \hspace*{4ex} {\tt common /addynvars\_diffkr/} &
977 \hspace*{4ex} is related to {\it DYNVARS.h} \\
978 \hspace*{4ex} {\tt common /addynvars\_kapgm/} &
979 \hspace*{4ex} is related to {\it DYNVARS.h} \\
980 \hspace*{4ex} {\tt common /adtr1\_r/} &
981 \hspace*{4ex} is related to {\it TR1.h} \\
982 \hspace*{4ex} {\tt common /adffields/} &
983 \hspace*{4ex} is related to {\it FFIELDS.h}\\
984 \end{tabular}
985
986 Note that if the structure of the common block changes in the
987 above header files of the forward code, the structure
988 of the adjoint common blocks will change accordingly.
989 Thus, it has to be made sure that the structure of the
990 adjoint common block in the hand-written file {\it adcommon.h}
991 complies with the automatically generated adjoint common blocks
992 in {\it adjoint\_model.F}.
993
994 \subsubsection{File {\it tamc.h}}
995
996 This routine contains the dimensions for TAMC checkpointing.
997 %
998 \begin{itemize}
999 %
1000 \item {\tt \#ifdef ALLOW\_TAMC\_CHECKPOINTING} \\
1001 3-level checkpointing is enabled, i.e. the timestepping
1002 is divided into three different levels (see Section \ref{???}).
1003 The model state of the outermost ({\tt nchklev\_3}) and the
1004 intermediate ({\tt nchklev\_2}) timestepping loop are stored to file
1005 (handled in {\it the\_main\_loop}).
1006 The innermost loop ({\tt nchklev\_1})
1007 avoids I/O by storing all required variables
1008 to common blocks. This storing may also be necessary if
1009 no checkpointing is chosen
1010 (nonlinear functions, if-statements, iterative loops, ...).
1011 In the present example the dimensions are chosen as follows: \\
1012 \hspace*{4ex} {\tt nchklev\_1 = 36 } \\
1013 \hspace*{4ex} {\tt nchklev\_2 = 30 } \\
1014 \hspace*{4ex} {\tt nchklev\_3 = 60 } \\
1015 To guarantee that the checkpointing intervals span the entire
1016 integration period the following relation must be satisfied: \\
1017 \hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\
1018 where {\tt nTimeSteps} is either specified in {\it data}
1019 or computed via \\
1020 \hspace*{4ex} {\tt nTimeSteps = (endTime-startTime)/deltaTClock }.
1021 %
1022 \item {\tt \#undef ALLOW\_TAMC\_CHECKPOINTING} \\
1023 No checkpointing is enabled.
1024 In this case the relevant counter is {\tt nchklev\_0}.
1025 Similar to above, the following relation has to be satisfied \\
1026 \hspace*{4ex} {\tt nchklev\_0 $ \ge $ nTimeSteps}.
1027 %
1028 \end{itemize}
1029
1030 The following parameters may be worth describing: \\
1031 %
1032 \hspace*{4ex} {\tt isbyte} \\
1033 \hspace*{4ex} {\tt maxpass} \\
1034 ~
1035
1036 \subsubsection{File {\it makefile}}
1037
1038 This file contains all relevant paramter flags and
1039 lists to run TAMC or TAF.
1040 It is assumed that TAMC is available to you, either locally,
1041 being installed on your network, or remotely through the 'TAMC Utility'.
1042 TAMC is called with the command {\tt tamc} followed by a
1043 number of options. They are described in detail in the
1044 TAMC manual \cite{gie:99}.
1045 Here we briefly discuss the main flags used in the {\it makefile}
1046 %
1047 \begin{itemize}
1048 \item [{\tt tamc}] {\tt
1049 -input <variable names>
1050 -output <variable name> -r4 ... \\
1051 -toplevel <S/R name> -reverse <file names>
1052 }
1053 \end{itemize}
1054 %
1055 \begin{itemize}
1056 %
1057 \item {\tt -toplevel <S/R name>} \\
1058 Name of the toplevel routine, with respect to which the
1059 control flow analysis is performed.
1060 %
1061 \item {\tt -input <variable names>} \\
1062 List of independent variables $ u $ with respect to which the
1063 dependent variable $ J $ is differentiated.
1064 %
1065 \item {\tt -output <variable name>} \\
1066 Dependent variable $ J $ which is to be differentiated.
1067 %
1068 \item {\tt -reverse <file names>} \\
1069 Adjoint code is generated to compute the sensitivity of an
1070 independent variable w.r.t. many dependent variables.
1071 In the discussion of Section ???
1072 the generated adjoint top-level routine computes the product
1073 of the transposed Jacobian matrix $ M^T $ times
1074 the gradient vector $ \nabla_v J $.
1075 \\
1076 {\tt <file names>} refers to the list of files {\it .f} which are to be
1077 analyzed by TAMC. This list is generally smaller than the full list
1078 of code to be compiled. The files not contained are either
1079 above the top-level routine (some initialisations), or are
1080 deliberately hidden from TAMC, either because hand-written
1081 adjoint routines exist, or the routines must not (or don't have to)
1082 be differentiated. For each routine which is part of the flow tree
1083 of the top-level routine, but deliberately hidden from TAMC
1084 (or for each package which contains such routines),
1085 a corresponding file {\it .flow} exists containing flow directives
1086 for TAMC.
1087 %
1088 \item {\tt -r4} \\
1089 ~
1090 %
1091 \end{itemize}
1092
1093
1094 \subsubsection{File {\it data}}
1095
1096 \subsubsection{File {\it data.cost}}
1097
1098 \subsubsection{File {\it data.ctrl}}
1099
1100 \subsubsection{File {\it data.gmredi}}
1101
1102 \subsubsection{File {\it data.grdchk}}
1103
1104 \subsubsection{File {\it data.optim}}
1105
1106 \subsubsection{File {\it data.pkg}}
1107
1108 \subsubsection{File {\it eedata}}
1109
1110 \subsubsection{File {\it topog.bin}}
1111
1112 \subsubsection{File {\it windx.bin, windy.bin}}
1113
1114 \subsubsection{File {\it salt.bin, theta.bin}}
1115
1116 \subsubsection{File {\it SSS.bin, SST.bin}}
1117
1118 \subsubsection{File {\it pickup*}}
1119
1120 \subsection{Compiling the model and its adjoint}
1121
1122 \newpage
1123
1124 %**********************************************************************
1125 \section{TLM and ADM generation in general}
1126 \label{sec_ad_setup_gen}
1127 %**********************************************************************
1128
1129 In this section we describe in a general fashion
1130 the parts of the code that are relevant for automatic
1131 differentiation using the software tool TAMC.
1132
1133 \begin{figure}[b!]
1134 \input{part5/doc_ad_the_model}
1135 \caption{~}
1136 \label{fig:adthemodel}
1137 \end{figure}
1138
1139 The basic flow is depicted in \reffig{adthemodel}.
1140 If the option {\tt ALLOW\_AUTODIFF\_TAMC} is defined, the driver routine
1141 {\it the\_model\_main}, instead of calling {\it the\_main\_loop},
1142 invokes the adjoint of this routine, {\it adthe\_main\_loop},
1143 which is the toplevel routine in terms of reverse mode computation.
1144 The routine {\it adthe\_main\_loop} has been generated using TAMC.
1145 It contains both the forward integration of the full model,
1146 any additional storing that is required for efficient checkpointing,
1147 and the reverse integration of the adjoint model.
1148 The structure of {\it adthe\_main\_loop} has been strongly
1149 simplified for clarification; in particular, no checkpointing
1150 procedures are shown here.
1151 Prior to the call of {\it adthe\_main\_loop}, the routine
1152 {\it ctrl\_unpack} is invoked to unpack the control vector,
1153 and following that call, the routine {\it ctrl\_pack}
1154 is invoked to pack the control vector
1155 (cf. Section \ref{section_ctrl}).
1156 If gradient checks are to be performed, the option
1157 {\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case
1158 the driver routine {\it grdchk\_main} is called after
1159 the gradient has been computed via the adjoint
1160 (cf. Section \ref{section_grdchk}).
1161
1162 \subsection{The cost function (dependent variable)
1163 \label{section_cost}}
1164
1165 The cost function $ {\cal J} $ is referred to as the {\sf dependent variable}.
1166 It is a function of the input variables $ \vec{u} $ via the composition
1167 $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.
1168 The input is referred to as the
1169 {\sf independent variables} or {\sf control variables}.
1170 All aspects relevant to the treatment of the cost function $ {\cal J} $
1171 (parameter setting, initialisation, accumulation,
1172 final evaluation), are controlled by the package {\it pkg/cost}.
1173
1174 \begin{figure}[h!]
1175 \input{part5/doc_cost_flow}
1176 \caption{~}
1177 \label{fig:costflow}
1178 \end{figure}
1179
1180 \subsubsection{genmake and CPP options}
1181 %
1182 \begin{itemize}
1183 %
1184 \item
1185 \fbox{
1186 \begin{minipage}{12cm}
1187 {\it genmake}, {\it CPP\_OPTIONS.h}, {\it ECCO\_CPPOPTIONS.h}
1188 \end{minipage}
1189 }
1190 \end{itemize}
1191 %
1192 The directory {\it pkg/cost} can be included to the
1193 compile list in 3 different ways (cf. Section \ref{???}):
1194 %
1195 \begin{enumerate}
1196 %
1197 \item {\it genmake}: \\
1198 Change the default settings in the file {\it genmake} by adding
1199 {\bf cost} to the {\bf enable} list (not recommended).
1200 %
1201 \item {\it .genmakerc}: \\
1202 Customize the settings of {\bf enable}, {\bf disable} which are
1203 appropriate for your experiment in the file {\it .genmakerc}
1204 and add the file to your compile directory.
1205 %
1206 \item genmake-options: \\
1207 Call {\it genmake} with the option
1208 {\tt genmake -enable=cost}.
1209 %
1210 \end{enumerate}
1211 The basic CPP option to enable the cost function is {\bf ALLOW\_COST}.
1212 Each specific cost function contribution has its own option.
1213 For the present example the option is {\bf ALLOW\_COST\_TRACER}.
1214 All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}
1215 Since the cost function is usually used in conjunction with
1216 automatic differentiation, the CPP option
1217 {\bf ALLOW\_ADJOINT\_RUN} should be defined
1218 (file {\it CPP\_OPTIONS.h}).
1219
1220 \subsubsection{Initialisation}
1221 %
1222 The initialisation of the {\it cost} package is readily enabled
1223 as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined.
1224 %
1225 \begin{itemize}
1226 %
1227 \item
1228 \fbox{
1229 \begin{minipage}{12cm}
1230 Parameters: {\it cost\_readparms}
1231 \end{minipage}
1232 }
1233 \\
1234 This S/R
1235 reads runtime flags and parameters from file {\it data.cost}.
1236 For the present example the only relevant parameter read
1237 is {\bf mult\_tracer}. This multiplier enables different
1238 cost function contributions to be switched on
1239 ( = 1.) or off ( = 0.) at runtime.
1240 For more complex cost functions which involve model vs. data
1241 misfits, the corresponding data filenames and data
1242 specifications (start date and time, period, ...) are read
1243 in this S/R.
1244 %
1245 \item
1246 \fbox{
1247 \begin{minipage}{12cm}
1248 Variables: {\it cost\_init}
1249 \end{minipage}
1250 }
1251 \\
1252 This S/R
1253 initialises the different cost function contributions.
1254 The contribtion for the present example is {\bf objf\_tracer}
1255 which is defined on each tile (bi,bj).
1256 %
1257 \end{itemize}
1258 %
1259 \subsubsection{Accumulation}
1260 %
1261 \begin{itemize}
1262 %
1263 \item
1264 \fbox{
1265 \begin{minipage}{12cm}
1266 {\it cost\_tile}, {\it cost\_tracer}
1267 \end{minipage}
1268 }
1269 \end{itemize}
1270 %
1271 The 'driver' routine
1272 {\it cost\_tile} is called at the end of each time step.
1273 Within this 'driver' routine, S/R are called for each of
1274 the chosen cost function contributions.
1275 In the present example ({\bf ALLOW\_COST\_TRACER}),
1276 S/R {\it cost\_tracer} is called.
1277 It accumulates {\bf objf\_tracer} according to eqn. (\ref{???}).
1278 %
1279 \subsubsection{Finalize all contributions}
1280 %
1281 \begin{itemize}
1282 %
1283 \item
1284 \fbox{
1285 \begin{minipage}{12cm}
1286 {\it cost\_final}
1287 \end{minipage}
1288 }
1289 \end{itemize}
1290 %
1291 At the end of the forward integration S/R {\it cost\_final}
1292 is called. It accumulates the total cost function {\bf fc}
1293 from each contribution and sums over all tiles:
1294 \begin{equation}
1295 {\cal J} \, = \,
1296 {\rm fc} \, = \,
1297 {\rm mult\_tracer} \sum_{bi,\,bj}^{nSx,\,nSy}
1298 {\rm objf\_tracer}(bi,bj) \, + \, ...
1299 \end{equation}
1300 %
1301 The total cost function {\bf fc} will be the
1302 'dependent' variable in the argument list for TAMC, i.e.
1303 \begin{verbatim}
1304 tamc -output 'fc' ...
1305 \end{verbatim}
1306
1307 %%%% \end{document}
1308
1309 \begin{figure}
1310 \input{part5/doc_ad_the_main}
1311 \caption{~}
1312 \label{fig:adthemain}
1313 \end{figure}
1314
1315 \subsection{The control variables (independent variables)
1316 \label{section_ctrl}}
1317
1318 The control variables are a subset of the model input
1319 (initial conditions, boundary conditions, model parameters).
1320 Here we identify them with the variable $ \vec{u} $.
1321 All intermediate variables whose derivative w.r.t. control
1322 variables do not vanish are called {\sf active variables}.
1323 All subroutines whose derivative w.r.t. the control variables
1324 don't vanish are called {\sf active routines}.
1325 Read and write operations from and to file can be viewed
1326 as variable assignments. Therefore, files to which
1327 active variables are written and from which active variables
1328 are read are called {\sf active files}.
1329 All aspects relevant to the treatment of the control variables
1330 (parameter setting, initialisation, perturbation)
1331 are controled by the package {\it pkg/ctrl}.
1332
1333 \begin{figure}[h!]
1334 \input{part5/doc_ctrl_flow}
1335 \caption{~}
1336 \label{fig:ctrlflow}
1337 \end{figure}
1338
1339 \subsubsection{genmake and CPP options}
1340 %
1341 \begin{itemize}
1342 %
1343 \item
1344 \fbox{
1345 \begin{minipage}{12cm}
1346 {\it genmake}, {\it CPP\_OPTIONS.h}, {\it ECCO\_CPPOPTIONS.h}
1347 \end{minipage}
1348 }
1349 \end{itemize}
1350 %
1351 To enable the directory to be included to the compile list,
1352 {\bf ctrl} has to be added to the {\bf enable} list in
1353 {\it .genmakerc} (or {\it genmake} itself).
1354 Each control variable is enabled via its own CPP option
1355 in {\it ECCO\_CPPOPTIONS.h}.
1356
1357 \subsubsection{Initialisation}
1358 %
1359 \begin{itemize}
1360 %
1361 \item
1362 \fbox{
1363 \begin{minipage}{12cm}
1364 Parameters: {\it ctrl\_readparms}
1365 \end{minipage}
1366 }
1367 \\
1368 %
1369 This S/R
1370 reads runtime flags and parameters from file {\it data.ctrl}.
1371 For the present example the file contains the file names
1372 of each control variable that is used.
1373 In addition, the number of wet points for each control
1374 variable and the net dimension of the space of control
1375 variables (counting wet points only) {\bf nvarlength}
1376 is determined.
1377 Masks for wet points for each tile {\bf (bi,\,bj)}
1378 and vertical layer {\bf k} are generated for the three
1379 relevant categories on the C-grid:
1380 {\bf nWetCtile} for tracer fields,
1381 {\bf nWetWtile} for zonal velocity fields,
1382 {\bf nWetStile} for meridional velocity fields.
1383 %
1384 \item
1385 \fbox{
1386 \begin{minipage}{12cm}
1387 Control variables, control vector,
1388 and their gradients: {\it ctrl\_unpack}
1389 \end{minipage}
1390 }
1391 \\
1392 %
1393 Two important issues related to the handling of the control
1394 variables in the MITGCM need to be addressed.
1395 First, in order to save memory, the control variable arrays
1396 are not kept in memory, but rather read from file and added
1397 to the initial fields during the model initialisation phase.
1398 Similarly, the corresponding adjoint fields which represent
1399 the gradient of the cost function w.r.t. the control variables
1400 are written to file at the end of the adjoint integration.
1401 Second, in addition to the files holding the 2-dim. and 3-dim.
1402 control variables and the corresponding cost gradients,
1403 a 1-dim. {\sf control vector}
1404 and {\sf gradient vector} are written to file. They contain
1405 only the wet points of the control variables and the corresponding
1406 gradient.
1407 This leads to a significant data compression.
1408 Furthermore, an option is available
1409 ({\tt ALLOW\_NONDIMENSIONAL\_CONTROL\_IO}) to
1410 non-dimensionalise the control and gradient vector,
1411 which otherwise would contain different pieces of different
1412 magnitudes and units.
1413 Finally, the control and gradient vector can be passed to a
1414 minimization routine if an update of the control variables
1415 is sought as part of a minimization exercise.
1416
1417 The files holding fields and vectors of the control variables
1418 and gradient are generated and initialised in S/R {\it ctrl\_unpack}.
1419 %
1420 \end{itemize}
1421
1422 \subsubsection{Perturbation of the independent variables}
1423 %
1424 The dependency flow for differentiation w.r.t. the controls
1425 starts with adding a perturbation onto the input variable,
1426 thus defining the independent or control variables for TAMC.
1427 Three types of controls may be considered:
1428 %
1429 \begin{itemize}
1430 %
1431 \item
1432 \fbox{
1433 \begin{minipage}{12cm}
1434 {\it ctrl\_map\_ini} (initial value sensitivity):
1435 \end{minipage}
1436 }
1437 \\
1438 %
1439 Consider as an example the initial tracer distribution
1440 {\bf tr1} as control variable.
1441 After {\bf tr1} has been initialised in
1442 {\it ini\_tr1} (dynamical variables such as
1443 temperature and salinity are initialised in {\it ini\_fields}),
1444 a perturbation anomaly is added to the field in S/R
1445 {\it ctrl\_map\_ini}
1446 %
1447 \begin{equation}
1448 \begin{split}
1449 u & = \, u_{[0]} \, + \, \Delta u \\
1450 {\bf tr1}(...) & = \, {\bf tr1_{ini}}(...) \, + \, {\bf xx\_tr1}(...)
1451 \label{perturb}
1452 \end{split}
1453 \end{equation}
1454 %
1455 {\bf xx\_tr1} is a 3-dim. global array
1456 holding the perturbation. In the case of a simple
1457 sensitivity study this array is identical to zero.
1458 However, it's specification is essential in the context
1459 of automatic differentiation since TAMC
1460 treats the corresponding line in the code symbolically
1461 when determining the differentiation chain and its origin.
1462 Thus, the variable names are part of the argument list
1463 when calling TAMC:
1464 %
1465 \begin{verbatim}
1466 tamc -input 'xx_tr1 ...' ...
1467 \end{verbatim}
1468 %
1469 Now, as mentioned above, the MITGCM avoids maintaining
1470 an array for each control variable by reading the
1471 perturbation to a temporary array from file.
1472 To ensure the symbolic link to be recognized by TAMC, a scalar
1473 dummy variable {\bf xx\_tr1\_dummy} is introduced
1474 and an 'active read' routine of the adjoint support
1475 package {\it pkg/autodiff} is invoked.
1476 The read-procedure is tagged with the variable
1477 {\bf xx\_tr1\_dummy} enabbling TAMC to recognize the
1478 initialisation of the perturbation.
1479 The modified call of TAMC thus reads
1480 %
1481 \begin{verbatim}
1482 tamc -input 'xx_tr1_dummy ...' ...
1483 \end{verbatim}
1484 %
1485 and the modified operation to (\ref{perturb})
1486 in the code takes on the form
1487 %
1488 \begin{verbatim}
1489 call active_read_xyz(
1490 & ..., tmpfld3d, ..., xx_tr1_dummy, ... )
1491
1492 tr1(...) = tr1(...) + tmpfld3d(...)
1493 \end{verbatim}
1494 %
1495 Note, that reading an active variable corresponds
1496 to a variable assignment. Its derivative corresponds
1497 to a write statement of the adjoint variable.
1498 The 'active file' routines have been designed
1499 to support active read and corresponding adjoint active write
1500 operations (and vice versa).
1501 %
1502 \item
1503 \fbox{
1504 \begin{minipage}{12cm}
1505 {\it ctrl\_map\_forcing} (boundary value sensitivity):
1506 \end{minipage}
1507 }
1508 \\
1509 %
1510 The handling of boundary values as control variables
1511 proceeds exactly analogous to the initial values
1512 with the symbolic perturbation taking place in S/R
1513 {\it ctrl\_map\_forcing}.
1514 Note however an important difference:
1515 Since the boundary values are time dependent with a new
1516 forcing field applied at each time steps,
1517 the general problem may be thought of as
1518 a new control variable at each time step
1519 (or, if the perturbation is averaged over a certain period,
1520 at each $ N $ timesteps), i.e.
1521 \[
1522 u_{\rm forcing} \, = \,
1523 \{ \, u_{\rm forcing} ( t_n ) \, \}_{
1524 n \, = \, 1, \ldots , {\rm nTimeSteps} }
1525 \]
1526 %
1527 In the current example an equilibrium state is considered,
1528 and only an initial perturbation to
1529 surface forcing is applied with respect to the
1530 equilibrium state.
1531 A time dependent treatment of the surface forcing is
1532 implemented in the ECCO environment, involving the
1533 calendar ({\it cal}~) and external forcing ({\it exf}~) packages.
1534 %
1535 \item
1536 \fbox{
1537 \begin{minipage}{12cm}
1538 {\it ctrl\_map\_params} (parameter sensitivity):
1539 \end{minipage}
1540 }
1541 \\
1542 %
1543 This routine is not yet implemented, but would proceed
1544 proceed along the same lines as the initial value sensitivity.
1545 The mixing parameters {\bf diffkr} and {\bf kapgm}
1546 are currently added as controls in {\it ctrl\_map\_ini.F}.
1547 %
1548 \end{itemize}
1549 %
1550
1551 \subsubsection{Output of adjoint variables and gradient}
1552 %
1553 Several ways exist to generate output of adjoint fields.
1554 %
1555 \begin{itemize}
1556 %
1557 \item
1558 \fbox{
1559 \begin{minipage}{12cm}
1560 {\it ctrl\_map\_ini, ctrl\_map\_forcing}:
1561 \end{minipage}
1562 }
1563 \\
1564 \begin{itemize}
1565 %
1566 \item {\bf xx\_...}: the control variable fields \\
1567 Before the forward integration, the control
1568 variables are read from file {\bf xx\_ ...} and added to
1569 the model field.
1570 %
1571 \item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient
1572 $ \nabla _{u}{\cal J} $ for each control variable \\
1573 After the adjoint integration the corresponding adjoint
1574 variables are written to {\bf adxx\_ ...}.
1575 %
1576 \end{itemize}
1577 %
1578 \item
1579 \fbox{
1580 \begin{minipage}{12cm}
1581 {\it ctrl\_unpack, ctrl\_pack}:
1582 \end{minipage}
1583 }
1584 \\
1585 %
1586 \begin{itemize}
1587 %
1588 \item {\bf vector\_ctrl}: the control vector \\
1589 At the very beginning of the model initialisation,
1590 the updated compressed control vector is read (or initialised)
1591 and distributed to 2-dim. and 3-dim. control variable fields.
1592 %
1593 \item {\bf vector\_grad}: the gradient vector \\
1594 At the very end of the adjoint integration,
1595 the 2-dim. and 3-dim. adjoint variables are read,
1596 compressed to a single vector and written to file.
1597 %
1598 \end{itemize}
1599 %
1600 \item
1601 \fbox{
1602 \begin{minipage}{12cm}
1603 {\it addummy\_in\_stepping}:
1604 \end{minipage}
1605 }
1606 \\
1607 In addition to writing the gradient at the end of the
1608 forward/adjoint integration, many more adjoint variables
1609 of the model state
1610 at intermediate times can be written using S/R
1611 {\it addummy\_in\_stepping}.
1612 This routine is part of the adjoint support package
1613 {\it pkg/autodiff} (cf.f. below).
1614 To be part of the adjoint code, the corresponding S/R
1615 {\it dummy\_in\_stepping} has to be called in the forward
1616 model (S/R {\it the\_main\_loop}) at the appropriate place.
1617
1618 {\it dummy\_in\_stepping} is essentially empty,
1619 the corresponding adjoint routine is hand-written rather
1620 than generated automatically.
1621 Appropriate flow directives ({\it dummy\_in\_stepping.flow})
1622 ensure that TAMC does not automatically
1623 generate {\it addummy\_in\_stepping} by trying to differentiate
1624 {\it dummy\_in\_stepping}, but instead refers to
1625 the hand-written routine.
1626
1627 {\it dummy\_in\_stepping} is called in the forward code
1628 at the beginning of each
1629 timestep, before the call to {\it dynamics}, thus ensuring
1630 that {\it addummy\_in\_stepping} is called at the end of
1631 each timestep in the adjoint calculation, after the call to
1632 {\it addynamics}.
1633
1634 {\it addummy\_in\_stepping} includes the header files
1635 {\it adcommon.h}.
1636 This header file is also hand-written. It contains
1637 the common blocks
1638 {\bf /addynvars\_r/}, {\bf /addynvars\_cd/},
1639 {\bf /addynvars\_diffkr/}, {\bf /addynvars\_kapgm/},
1640 {\bf /adtr1\_r/}, {\bf /adffields/},
1641 which have been extracted from the adjoint code to enable
1642 access to the adjoint variables.
1643 %
1644 \end{itemize}
1645
1646
1647 \subsubsection{Control variable handling for
1648 optimization applications}
1649
1650 In optimization mode the cost function $ {\cal J}(u) $ is sought
1651 to be minimized with respect to a set of control variables
1652 $ \delta {\cal J} \, = \, 0 $, in an iterative manner.
1653 The gradient $ \nabla _{u}{\cal J} |_{u_{[k]}} $ together
1654 with the value of the cost function itself $ {\cal J}(u_{[k]}) $
1655 at iteration step $ k $ serve
1656 as input to a minimization routine (e.g. quasi-Newton method,
1657 conjugate gradient, ... \cite{gil_lem:89})
1658 to compute an update in the
1659 control variable for iteration step $k+1$
1660 \[
1661 u_{[k+1]} \, = \, u_{[0]} \, + \, \Delta u_{[k+1]}
1662 \quad \mbox{satisfying} \quad
1663 {\cal J} \left( u_{[k+1]} \right) \, < \, {\cal J} \left( u_{[k]} \right)
1664 \]
1665 $ u_{[k+1]} $ then serves as input for a forward/adjoint run
1666 to determine $ {\cal J} $ and $ \nabla _{u}{\cal J} $ at iteration step
1667 $ k+1 $.
1668 Tab. \ref{???} sketches the flow between forward/adjoint model
1669 and the minimization routine.
1670
1671 \begin{eqnarray*}
1672 \scriptsize
1673 \begin{array}{ccccc}
1674 u_{[0]} \,\, , \,\, \Delta u_{[k]} & ~ & ~ & ~ & ~ \\
1675 {\Big\downarrow}
1676 & ~ & ~ & ~ & ~ \\
1677 ~ & ~ & ~ & ~ & ~ \\
1678 \hline
1679 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1680 \multicolumn{1}{|c}{
1681 u_{[k]} = u_{[0]} + \Delta u_{[k]}} &
1682 \stackrel{\bf forward}{\bf \longrightarrow} &
1683 v_{[k]} = M \left( u_{[k]} \right) &
1684 \stackrel{\bf forward}{\bf \longrightarrow} &
1685 \multicolumn{1}{c|}{
1686 {\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\
1687 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1688 \hline
1689 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1690 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{{\Big\downarrow}} \\
1691 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1692 \hline
1693 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1694 \multicolumn{1}{|c}{
1695 \nabla_u {\cal J}_{[k]} (\delta {\cal J}) =
1696 T^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} &
1697 \stackrel{\bf adjoint}{\mathbf \longleftarrow} &
1698 ad \, v_{[k]} (\delta {\cal J}) =
1699 \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) &
1700 \stackrel{\bf adjoint}{\mathbf \longleftarrow} &
1701 \multicolumn{1}{c|}{ ad \, {\cal J} = \delta {\cal J}} \\
1702 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1703 \hline
1704 ~ & ~ & ~ & ~ & ~ \\
1705 \hspace*{15ex}{\Bigg\downarrow}
1706 \quad {\cal J}_{[k]}, \quad \nabla_u {\cal J}_{[k]}
1707 & ~ & ~ & ~ & ~ \\
1708 ~ & ~ & ~ & ~ & ~ \\
1709 \hline
1710 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1711 \multicolumn{1}{|c}{
1712 {\cal J}_{[k]} \,\, , \,\, \nabla_u {\cal J}_{[k]}} &
1713 {\mathbf \longrightarrow} & \text{\bf minimisation} &
1714 {\mathbf \longrightarrow} &
1715 \multicolumn{1}{c|}{ \Delta u_{[k+1]}} \\
1716 \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1717 \hline
1718 ~ & ~ & ~ & ~ & ~ \\
1719 ~ & ~ & ~ & ~ & \Big\downarrow \\
1720 ~ & ~ & ~ & ~ & \Delta u_{[k+1]} \\
1721 \end{array}
1722 \end{eqnarray*}
1723
1724 The routines {\it ctrl\_unpack} and {\it ctrl\_pack} provide
1725 the link between the model and the minimization routine.
1726 As described in Section \ref{???}
1727 the {\it unpack} and {\it pack} routines read and write
1728 control and gradient {\it vectors} which are compressed
1729 to contain only wet points, in addition to the full
1730 2-dim. and 3-dim. fields.
1731 The corresponding I/O flow looks as follows:
1732
1733 \vspace*{0.5cm}
1734
1735 {\scriptsize
1736 \begin{tabular}{ccccc}
1737 {\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\
1738 {\big\downarrow} & ~ & ~ & ~ & ~ \\
1739 \cline{1-1}
1740 \multicolumn{1}{|c|}{\it ctrl\_unpack} & ~ & ~ & ~ & ~ \\
1741 \cline{1-1}
1742 {\big\downarrow} & ~ & ~ & ~ & ~ \\
1743 \cline{3-3}
1744 \multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ &
1745 \multicolumn{1}{|c|}{~} & ~ & ~ \\
1746 \multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} &
1747 $\stackrel{\mbox{read}}{\longrightarrow}$ &
1748 \multicolumn{1}{|c|}{forward integration} & ~ & ~ \\
1749 \multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~}
1750 & ~ & ~ \\
1751 \cline{3-3}
1752 ~ & ~ & $\downarrow$ & ~ & ~ \\
1753 \cline{3-3}
1754 ~ & ~ &
1755 \multicolumn{1}{|c|}{~} & ~ &
1756 \multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$} \\
1757 ~ & ~ & \multicolumn{1}{|c|}{adjoint integration} &
1758 $\stackrel{\mbox{write}}{\longrightarrow}$ &
1759 \multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\
1760 ~ & ~ & \multicolumn{1}{|c|}{~}
1761 & ~ & \multicolumn{1}{l}{\bf \vdots} \\
1762 \cline{3-3}
1763 ~ & ~ & ~ & ~ & {\big\downarrow} \\
1764 \cline{5-5}
1765 ~ & ~ & ~ & ~ & \multicolumn{1}{|c|}{\it ctrl\_pack} \\
1766 \cline{5-5}
1767 ~ & ~ & ~ & ~ & {\big\downarrow} \\
1768 ~ & ~ & ~ & ~ & {\bf vector\_grad\_$<$k$>$ } \\
1769 \end{tabular}
1770 }
1771
1772 \vspace*{0.5cm}
1773
1774
1775 {\it ctrl\_unpack} reads the updated control vector
1776 {\bf vector\_ctrl\_$<$k$>$}.
1777 It distributes the different control variables to
1778 2-dim. and 3-dim. files {\it xx\_...$<$k$>$}.
1779 At the start of the forward integration the control variables
1780 are read from {\it xx\_...$<$k$>$} and added to the
1781 field.
1782 Correspondingly, at the end of the adjoint integration
1783 the adjoint fields are written
1784 to {\it adxx\_...$<$k$>$}, again via the active file routines.
1785 Finally, {\it ctrl\_pack} collects all adjoint files
1786 and writes them to the compressed vector file
1787 {\bf vector\_grad\_$<$k$>$}.
1788
1789 \subsection{TLM and ADM generation via TAMC}
1790
1791
1792
1793 \subsection{Flow directives and adjoint support routines \label{section_flowdir}}
1794
1795 \subsection{Store directives and checkpointing \label{section_checkpointing}}
1796
1797 \subsection{Gradient checks \label{section_grdchk}}
1798
1799 \subsection{Second derivative generation via TAMC}
1800
1801 \section{Example of adjoint code}

  ViewVC Help
Powered by ViewVC 1.1.22