/[MITgcm]/MITgcm/pkg/grdchk/grdchk_main.F
ViewVC logotype

Contents of /MITgcm/pkg/grdchk/grdchk_main.F

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


Revision 1.42 - (show annotations) (download)
Sat Aug 25 14:18:51 2012 UTC (11 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64, checkpoint63s, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.41: +1 -3 lines
Always includes "g_cost.h": should not have conditional included files depending
on any option from AD_CONFIG.h (ALLOW_ADJOINT_RUN or ALLOW_TANGENTLINEAR_RUN):
this would produce wrong dependency in Makefile.

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

  ViewVC Help
Powered by ViewVC 1.1.22