/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_common/grdchk_main.F
ViewVC logotype

Contents of /MITgcm_contrib/heimbach/OpenAD/code_common/grdchk_main.F

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


Revision 1.3 - (show annotations) (download)
Wed Aug 22 22:35:45 2012 UTC (12 years, 11 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
FILE REMOVED
changes roled in main version

1
2 C $Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_common/grdchk_main.F,v 1.2 2011/09/26 20:22:18 utke Exp $
3 C $Name: $
4
5 #include "AD_CONFIG.h"
6 #include "CPP_OPTIONS.h"
7
8 CBOI
9 C
10 C !TITLE: GRADIENT CHECK
11 C !AUTHORS: mitgcm developers ( support@mitgcm.org )
12 C !AFFILIATION: Massachussetts Institute of Technology
13 C !DATE:
14 C !INTRODUCTION: gradient check package
15 c \bv
16 c Compare the gradients calculated by the adjoint model with
17 c finite difference approximations.
18 c
19 C !CALLING SEQUENCE:
20 c
21 c the_model_main
22 c |
23 c |-- ctrl_unpack
24 c |-- adthe_main_loop - unperturbed cost function and
25 c |-- ctrl_pack adjoint gradient are computed here
26 c |
27 c |-- grdchk_main
28 c |
29 c |-- grdchk_init
30 c |-- do icomp=... - loop over control vector elements
31 c |
32 c |-- grdchk_loc - determine location of icomp on grid
33 c |
34 c |-- grdchk_getxx - get control vector component from file
35 c | perturb it and write back to file
36 c |-- grdchk_getadxx - get gradient component calculated
37 c | via adjoint
38 c |-- the_main_loop - forward run and cost evaluation
39 c | with perturbed control vector element
40 c |-- calculate ratio of adj. vs. finite difference gradient
41 c |
42 c |-- grdchk_setxx - Reset control vector element
43 c |
44 c |-- grdchk_print - print results
45 c \ev
46 CEOI
47
48 CBOP
49 C !ROUTINE: grdchk_main
50 C !INTERFACE:
51 subroutine grdchk_main( mythid )
52
53 C !DESCRIPTION: \bv
54 c ==================================================================
55 c SUBROUTINE grdchk_main
56 c ==================================================================
57 c o Compare the gradients calculated by the adjoint model with
58 c finite difference approximations.
59 c started: Christian Eckert eckert@mit.edu 24-Feb-2000
60 c continued&finished: heimbach@mit.edu: 13-Jun-2001
61 c changed: mlosch@ocean.mit.edu: 09-May-2002
62 c - added centered difference vs. 1-sided difference option
63 c - improved output format for readability
64 c - added control variable hFacC
65 c heimbach@mit.edu 24-Feb-2003
66 c - added tangent linear gradient checks
67 c - fixes for multiproc. gradient checks
68 c - added more control variables
69 c
70 c ==================================================================
71 c SUBROUTINE grdchk_main
72 c ==================================================================
73 C \ev
74
75 C !USES:
76 implicit none
77
78 c == global variables ==
79 #include "SIZE.h"
80 #include "EEPARAMS.h"
81 #include "PARAMS.h"
82 #include "grdchk.h"
83 #include "cost.h"
84 #include "ctrl.h"
85 #ifdef ALLOW_TANGENTLINEAR_RUN
86 #include "g_cost.h"
87 #endif
88
89 C !INPUT/OUTPUT PARAMETERS:
90 c == routine arguments ==
91 integer mythid
92
93 #ifdef ALLOW_GRDCHK
94 C !LOCAL VARIABLES:
95 c == local variables ==
96 integer myiter
97 _RL mytime
98 integer bi, itlo, ithi
99 integer bj, jtlo, jthi
100 integer i, imin, imax
101 integer j, jmin, jmax
102 integer k
103
104 integer icomp
105 integer ichknum
106 integer icvrec
107 integer jtile
108 integer itile
109 integer layer
110 integer obcspos
111 integer itilepos
112 integer jtilepos
113 integer icglo
114 integer itest
115 integer ierr
116 integer ierr_grdchk
117 _RL gfd
118 _RL fcref
119 _RL fcpertplus, fcpertminus
120 _RL ratio_ad
121 _RL ratio_ftl
122 _RL xxmemo_ref
123 _RL xxmemo_pert
124 _RL adxxmemo
125 _RL ftlxxmemo
126 _RL localEps
127 _RL grdchk_epsfac
128
129 _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130 _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
131 _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
132
133 CHARACTER*(MAX_LEN_MBUF) msgBuf
134
135 c == end of interface ==
136 CEOP
137
138 c-- Set the loop ranges.
139 jtlo = mybylo(mythid)
140 jthi = mybyhi(mythid)
141 itlo = mybxlo(mythid)
142 ithi = mybxhi(mythid)
143 jmin = 1
144 jmax = sny
145 imin = 1
146 imax = snx
147
148 print *, 'ph-check entering grdchk_main '
149
150 c-- initialise variables
151 call grdchk_init( mythid )
152
153 c-- Compute the adjoint model gradients.
154 c-- Compute the unperturbed cost function.
155 cph Gradient via adjoint has already been computed,
156 cph and so has unperturbed cost function,
157 cph assuming all xx_ fields are initialised to zero.
158
159 ierr = 0
160 ierr_grdchk = 0
161 adxxmemo = 0.
162 ftlxxmemo = 0.
163 #ifdef ALLOW_ADMTLM
164 fcref = objf_state_final(idep,jdep,1,1,1)
165 #else
166 fcref = fc%v
167 #endif
168
169 print *, 'ph-check fcref = ', fcref
170
171 do bj = jtlo, jthi
172 do bi = itlo, ithi
173 do k = 1, nr
174 do j = jmin, jmax
175 do i = imin, imax
176 tmpplot1(i,j,k,bi,bj) = 0. _d 0
177 tmpplot2(i,j,k,bi,bj) = 0. _d 0
178 tmpplot3(i,j,k,bi,bj) = 0. _d 0
179 end do
180 end do
181 end do
182 end do
183 end do
184
185 if ( useCentralDiff ) then
186 grdchk_epsfac = 2. _d 0
187 else
188 grdchk_epsfac = 1. _d 0
189 end if
190
191 WRITE(standardmessageunit,'(A)')
192 & 'grad-res -------------------------------'
193 WRITE(standardmessageunit,'(2a)')
194 & ' grad-res proc # i j k bi bj iobc',
195 & ' fc ref fc + eps fc - eps'
196 #ifdef ALLOW_TANGENTLINEAR_RUN
197 WRITE(standardmessageunit,'(2a)')
198 & ' grad-res proc # i j k bi bj iobc',
199 & ' tlm grad fd grad 1 - fd/tlm'
200 #else
201 WRITE(standardmessageunit,'(2a)')
202 & ' grad-res proc # i j k bi bj iobc',
203 & ' adj grad fd grad 1 - fd/adj'
204 #endif
205
206 c-- Compute the finite difference approximations.
207 c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
208 c-- gradient checks.
209
210 if ( nbeg .EQ. 0 )
211 & call grdchk_get_position( mythid )
212
213 do icomp = nbeg, nend, nstep
214
215 ichknum = (icomp - nbeg)/nstep + 1
216 adxxmemo = 0.
217
218 cph(
219 cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
220 cph-print & nbeg, icomp, ichknum
221 cph)
222 if (ichknum .le. maxgrdchecks ) then
223
224 c-- Determine the location of icomp on the grid.
225 if ( myProcId .EQ. grdchkwhichproc ) then
226 call grdchk_loc( icomp, ichknum,
227 & icvrec, itile, jtile, layer, obcspos,
228 & itilepos, jtilepos, icglo, itest, ierr,
229 & mythid )
230 cph(
231 cph-print print *, 'ph-grd ----- back from loc -----',
232 cph-print & icvrec, itilepos, jtilepos, layer, obcspos
233 cph)
234 endif
235 _BARRIER
236
237 c******************************************************
238 c-- (A): get gradient component calculated via adjoint
239 c******************************************************
240
241 c-- get gradient component calculated via adjoint
242 if ( myProcId .EQ. grdchkwhichproc .AND.
243 & ierr .EQ. 0 ) then
244 call grdchk_getadxx( icvrec,
245 & itile, jtile, layer,
246 & itilepos, jtilepos,
247 & adxxmemo, mythid )
248 endif
249 C-- Add a global-sum call so that all proc will get the adjoint gradient
250 _GLOBAL_SUM_RL( adxxmemo, myThid )
251 c _BARRIER
252
253 #ifdef ALLOW_TANGENTLINEAR_RUN
254 c******************************************************
255 c-- (B): Get gradient component g_fc from tangent linear run:
256 c******************************************************
257 c--
258 c-- 1. perturb control vector component: xx(i)=1.
259
260 if ( myProcId .EQ. grdchkwhichproc .AND.
261 & ierr .EQ. 0 ) then
262 localEps = 1. _d 0
263 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
264 & itile, jtile, layer,
265 & itilepos, jtilepos,
266 & xxmemo_ref, xxmemo_pert, localEps,
267 & mythid )
268 endif
269 _BARRIER
270
271 c--
272 c-- 2. perform tangent linear run
273 mytime = starttime
274 myiter = niter0
275 #ifdef ALLOW_ADMTLM
276 do k=1,4*Nr+1
277 do j=1,sny
278 do i=1,snx
279 g_objf_state_final(i,j,1,1,k) = 0.
280 enddo
281 enddo
282 enddo
283 #else
284 g_fc = 0.
285 #endif
286
287 call g_the_main_loop( mytime, myiter, mythid )
288 _BARRIER
289 #ifdef ALLOW_ADMTLM
290 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
291 #else
292 ftlxxmemo = g_fc
293 #endif
294
295 c--
296 c-- 3. reset control vector
297 if ( myProcId .EQ. grdchkwhichproc .AND.
298 & ierr .EQ. 0 ) then
299 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
300 & itile, jtile, layer,
301 & itilepos, jtilepos,
302 & xxmemo_ref, mythid )
303 endif
304 _BARRIER
305
306 #endif /* ALLOW_TANGENTLINEAR_RUN */
307
308
309 c******************************************************
310 c-- (C): Get gradient via finite difference perturbation
311 c******************************************************
312
313 c-- get control vector component from file
314 c-- perturb it and write back to file
315 c-- positive perturbation
316 localEps = abs(grdchk_eps)
317 if ( myProcId .EQ. grdchkwhichproc .AND.
318 & ierr .EQ. 0 ) then
319 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
320 & itile, jtile, layer,
321 & itilepos, jtilepos,
322 & xxmemo_ref, xxmemo_pert, localEps,
323 & mythid )
324 endif
325 _BARRIER
326
327 c-- forward run with perturbed control vector
328 mytime = starttime
329 myiter = niter0
330 call OpenAD_the_main_loop( mytime, myiter, mythid )
331 #ifdef ALLOW_ADMTLM
332 fcpertplus = objf_state_final(idep,jdep,1,1,1)
333 #else
334 fcpertplus = fc%v
335 #endif
336 print *, 'ph-check fcpertplus = ', fcpertplus
337 _BARRIER
338
339 c-- Reset control vector.
340 if ( myProcId .EQ. grdchkwhichproc .AND.
341 & ierr .EQ. 0 ) then
342 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
343 & itile, jtile, layer,
344 & itilepos, jtilepos,
345 & xxmemo_ref, mythid )
346 endif
347 _BARRIER
348
349 fcpertminus = fcref
350 print *, 'ph-check fcpertminus = ', fcpertminus
351
352 if ( useCentralDiff ) then
353
354 c-- get control vector component from file
355 c-- perturb it and write back to file
356 c-- repeat the proceedure for a negative perturbation
357 if ( myProcId .EQ. grdchkwhichproc .AND.
358 & ierr .EQ. 0 ) then
359 localEps = - abs(grdchk_eps)
360 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
361 & itile, jtile, layer,
362 & itilepos, jtilepos,
363 & xxmemo_ref, xxmemo_pert, localEps,
364 & mythid )
365 endif
366 _BARRIER
367
368 c-- forward run with perturbed control vector
369 mytime = starttime
370 myiter = niter0
371 call OpenAD_the_main_loop( mytime, myiter, mythid )
372 _BARRIER
373 #ifdef ALLOW_ADMTLM
374 fcpertminus = objf_state_final(idep,jdep,1,1,1)
375 #else
376 fcpertminus = fc%v
377 #endif
378
379 c-- Reset control vector.
380 if ( myProcId .EQ. grdchkwhichproc .AND.
381 & ierr .EQ. 0 ) then
382 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
383 & itile, jtile, layer,
384 & itilepos, jtilepos,
385 & xxmemo_ref, mythid )
386 endif
387 _BARRIER
388
389 c-- end of if useCentralDiff ...
390 end if
391
392 c******************************************************
393 c-- (D): calculate relative differences between gradients
394 c******************************************************
395
396 if ( grdchk_eps .eq. 0. ) then
397 gfd = (fcpertplus-fcpertminus)
398 else
399 gfd = (fcpertplus-fcpertminus)
400 & /(grdchk_epsfac*grdchk_eps)
401 endif
402
403 if ( adxxmemo .eq. 0. ) then
404 ratio_ad = abs( adxxmemo - gfd )
405 else
406 ratio_ad = 1. - gfd/adxxmemo
407 endif
408
409 if ( ftlxxmemo .eq. 0. ) then
410 ratio_ftl = abs( ftlxxmemo - gfd )
411 else
412 ratio_ftl = 1. - gfd/ftlxxmemo
413 endif
414
415 if ( myProcId .EQ. grdchkwhichproc .AND.
416 & ierr .EQ. 0 ) then
417
418 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
419 & = gfd
420 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
421 & = ratio_ad
422 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
423 & = ratio_ftl
424
425 fcrmem ( ichknum ) = fcref
426 fcppmem ( ichknum ) = fcpertplus
427 fcpmmem ( ichknum ) = fcpertminus
428 xxmemref ( ichknum ) = xxmemo_ref
429 xxmempert ( ichknum ) = xxmemo_pert
430 gfdmem ( ichknum ) = gfd
431 adxxmem ( ichknum ) = adxxmemo
432 ftlxxmem ( ichknum ) = ftlxxmemo
433 ratioadmem ( ichknum ) = ratio_ad
434 ratioftlmem ( ichknum ) = ratio_ftl
435
436 irecmem ( ichknum ) = icvrec
437 bimem ( ichknum ) = itile
438 bjmem ( ichknum ) = jtile
439 ilocmem ( ichknum ) = itilepos
440 jlocmem ( ichknum ) = jtilepos
441 klocmem ( ichknum ) = layer
442 iobcsmem ( ichknum ) = obcspos
443 icompmem ( ichknum ) = icomp
444 ichkmem ( ichknum ) = ichknum
445 itestmem ( ichknum ) = itest
446 ierrmem ( ichknum ) = ierr
447 icglomem ( ichknum ) = icglo
448
449 WRITE(standardmessageunit,'(A)')
450 & 'grad-res -------------------------------'
451 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
452 & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
453 & layer,itile,jtile,obcspos,
454 & fcref, fcpertplus, fcpertminus
455 #ifdef ALLOW_TANGENTLINEAR_RUN
456 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
457 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
458 & icompmem(ichknum),itestmem(ichknum),
459 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
460 & ftlxxmemo, gfd, ratio_ftl
461 #else
462 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
463 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
464 & icompmem(ichknum),itestmem(ichknum),
465 & bimem(ichknum),bjmem(ichknum),obcspos,
466 & adxxmemo, gfd, ratio_ad
467 #endif
468 endif
469 #ifdef ALLOW_TANGENTLINEAR_RUN
470 WRITE(msgBuf,'(A34,1PE24.14)')
471 & ' TLM precision_derivative_cost =', fcref
472 CALL PRINT_MESSAGE
473 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
474 WRITE(msgBuf,'(A34,1PE24.14)')
475 & ' TLM precision_derivative_grad =', ftlxxmemo
476 CALL PRINT_MESSAGE
477 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
478 #else
479 WRITE(msgBuf,'(A30,1PE22.14)')
480 & ' ADM ref_cost_function =', fcref
481 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
482 & SQUEEZE_RIGHT, myThid )
483 WRITE(msgBuf,'(A30,1PE22.14)')
484 & ' ADM adjoint_gradient =', adxxmemo
485 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
486 & SQUEEZE_RIGHT, myThid )
487 WRITE(msgBuf,'(A30,1PE22.14)')
488 & ' ADM finite-diff_grad =', gfd
489 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
490 & SQUEEZE_RIGHT, myThid )
491 #endif
492
493 print *, 'ph-grd ierr ---------------------------'
494 print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
495 & ', ichknum = ', ichknum
496
497 _BARRIER
498
499 c-- else of if ( ichknum ...
500 else
501 ierr_grdchk = -1
502
503 c-- end of if ( ichknum ...
504 endif
505
506 c-- end of do icomp = ...
507 enddo
508
509 if ( myProcId .EQ. grdchkwhichproc ) then
510 CALL WRITE_REC_XYZ_RL(
511 & 'grd_findiff' , tmpplot1, 1, 0, myThid)
512 CALL WRITE_REC_XYZ_RL(
513 & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
514 CALL WRITE_REC_XYZ_RL(
515 & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
516 endif
517
518 c-- Everyone has to wait for the component to be reset.
519 _BARRIER
520
521 c-- Print the results of the gradient check.
522 call grdchk_print( ichknum, ierr_grdchk, mythid )
523
524 #endif /* ALLOW_GRDCHK */
525
526 return
527 end

  ViewVC Help
Powered by ViewVC 1.1.22