/[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.41 - (show annotations) (download)
Fri Aug 24 18:43:41 2012 UTC (11 years, 8 months ago) by jmc
Branch: MAIN
Changes since 1.40: +360 -372 lines
- add a call to TURNOFF_MODEL_IO for Tangent-Linear run (in this case, call
  to this routine from cost_final has been dropped in g_cost_final)
- improve printed information (more explicit msg, no longer using "print *,"
  fix fcpertminus printed value).

1 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.40 2012/08/22 19:34:51 heimbach 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 #ifdef ALLOW_TANGENTLINEAR_RUN
84 #include "g_cost.h"
85 #endif
86
87 C !INPUT/OUTPUT PARAMETERS:
88 C == routine arguments ==
89 INTEGER myThid
90
91 #ifdef ALLOW_GRDCHK
92 C !LOCAL VARIABLES:
93 C == local variables ==
94 INTEGER myIter
95 _RL myTime
96 INTEGER bi, bj
97 INTEGER i, j, k
98 INTEGER iMin, iMax, jMin, jMax
99 PARAMETER( iMin = 1 , iMax = sNx , jMin = 1 , jMax = sNy )
100 INTEGER ioUnit
101 CHARACTER*(MAX_LEN_MBUF) msgBuf
102
103 INTEGER icomp
104 INTEGER ichknum
105 INTEGER icvrec
106 INTEGER jtile
107 INTEGER itile
108 INTEGER layer
109 INTEGER obcspos
110 INTEGER itilepos
111 INTEGER jtilepos
112 INTEGER icglo
113 INTEGER itest
114 INTEGER ierr
115 INTEGER ierr_grdchk
116 _RL gfd
117 _RL fcref
118 _RL fcpertplus, fcpertminus
119 _RL ratio_ad
120 _RL ratio_ftl
121 _RL xxmemo_ref
122 _RL xxmemo_pert
123 _RL adxxmemo
124 _RL ftlxxmemo
125 _RL localEps
126 _RL grdchk_epsfac
127 _RL tmpplot1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
128 _RL tmpplot2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
129 _RL tmpplot3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
130 C == end of interface ==
131 CEOP
132
133 ioUnit = standardMessageUnit
134 WRITE(msgBuf,'(A)')
135 &'// ======================================================='
136 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
137 WRITE(msgBuf,'(A)') '// Gradient-check starts (grdchk_main)'
138 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
139 WRITE(msgBuf,'(A)')
140 &'// ======================================================='
141 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
142
143 #ifdef ALLOW_TANGENTLINEAR_RUN
144 C-- prevent writing output multiple times
145 C note: already called in AD run so that only needed for TLM run
146 CALL TURNOFF_MODEL_IO( 0, myThid )
147 #endif
148
149 C-- initialise variables
150 CALL GRDCHK_INIT( myThid )
151
152 C-- Compute the adjoint model gradients.
153 C-- Compute the unperturbed cost function.
154 Cph Gradient via adjoint has already been computed,
155 Cph and so has unperturbed cost function,
156 Cph assuming all xx_ fields are initialised to zero.
157
158 ierr = 0
159 ierr_grdchk = 0
160 adxxmemo = 0.
161 ftlxxmemo = 0.
162 #if (defined (ALLOW_ADMTLM))
163 fcref = objf_state_final(idep,jdep,1,1,1)
164 #elif (defined (ALLOW_AUTODIFF_OPENAD))
165 fcref = fc%v
166 #else
167 fcref = fc
168 #endif
169 WRITE(msgBuf,'(A,1PE22.14)')
170 & 'grdchk reference fc: fcref =', fcref
171 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
172
173 DO bj = myByLo(myThid), myByHi(myThid)
174 DO bi = myBxLo(myThid), myBxHi(myThid)
175 DO k = 1, Nr
176 DO j = jMin, jMax
177 DO i = iMin, iMax
178 tmpplot1(i,j,k,bi,bj) = 0. _d 0
179 tmpplot2(i,j,k,bi,bj) = 0. _d 0
180 tmpplot3(i,j,k,bi,bj) = 0. _d 0
181 ENDDO
182 ENDDO
183 ENDDO
184 ENDDO
185 ENDDO
186
187 IF ( useCentralDiff ) THEN
188 grdchk_epsfac = 2. _d 0
189 ELSE
190 grdchk_epsfac = 1. _d 0
191 ENDIF
192
193 WRITE(standardmessageunit,'(A)')
194 & 'grad-res -------------------------------'
195 WRITE(standardmessageunit,'(2a)')
196 & ' grad-res proc # i j k bi bj iobc',
197 & ' fc ref fc + eps fc - eps'
198 #ifdef ALLOW_TANGENTLINEAR_RUN
199 WRITE(standardmessageunit,'(2a)')
200 & ' grad-res proc # i j k bi bj iobc',
201 & ' tlm grad fd grad 1 - fd/tlm'
202 #else
203 WRITE(standardmessageunit,'(2a)')
204 & ' grad-res proc # i j k bi bj iobc',
205 & ' adj grad fd grad 1 - fd/adj'
206 #endif
207
208 C-- Compute the finite difference approximations.
209 C-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
210 C-- gradient checks.
211
212 IF ( nbeg .EQ. 0 ) CALL GRDCHK_GET_POSITION( myThid )
213
214 DO icomp = nbeg, nend, nstep
215
216 adxxmemo = 0.
217 ichknum = (icomp - nbeg)/nstep + 1
218
219 IF ( ichknum .LE. maxgrdchecks ) THEN
220 WRITE(msgBuf,'(A,I4,A)')
221 & '====== Starts gradient-check number', ichknum,
222 & ' (=ichknum) ======='
223 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
224
225 C-- Determine the location of icomp on the grid.
226 IF ( myProcId .EQ. grdchkwhichproc ) THEN
227 CALL grdchk_loc( icomp, ichknum,
228 & icvrec, itile, jtile, layer, obcspos,
229 & itilepos, jtilepos, icglo, itest, ierr,
230 & myThid )
231 ELSE
232 icvrec = 0
233 itile = 0
234 jtile = 0
235 layer = 0
236 obcspos = 0
237 itilepos = 0
238 jtilepos = 0
239 icglo = 0
240 itest = 0
241 ENDIF
242
243 C make sure that all procs have correct file records, so that useSingleCpuIO works
244 CALL GLOBAL_SUM_INT( icvrec , myThid )
245 CALL GLOBAL_SUM_INT( layer , myThid )
246 CALL GLOBAL_SUM_INT( ierr , myThid )
247
248 WRITE(msgBuf,'(A,3I5,A,2I4,A,I3,A,I4)')
249 & 'grdchk pos: i,j,k=', itilepos, jtilepos, layer,
250 & ' ; bi,bj=', itile, jtile, ' ; iobc=', obcspos,
251 & ' ; rec=', icvrec
252 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
253
254 C******************************************************
255 C-- (A): get gradient component calculated via adjoint
256 C******************************************************
257
258 C-- get gradient component calculated via adjoint
259 CALL GRDCHK_GETADXX( icvrec,
260 & itile, jtile, layer,
261 & itilepos, jtilepos,
262 & adxxmemo, ierr, myThid )
263 C-- Add a global-sum call so that all proc will get the adjoint gradient
264 _GLOBAL_SUM_RL( adxxmemo, myThid )
265
266 #ifdef ALLOW_TANGENTLINEAR_RUN
267 C******************************************************
268 C-- (B): Get gradient component g_fc from tangent linear run:
269 C******************************************************
270
271 C-- 1. perturb control vector component: xx(i)=1.
272
273 localEps = 1. _d 0
274 CALL GRDCHK_GETXX( icvrec, TANGENT_SIMULATION,
275 & itile, jtile, layer,
276 & itilepos, jtilepos,
277 & xxmemo_ref, xxmemo_pert, localEps,
278 & ierr, myThid )
279
280 C-- 2. perform tangent linear run
281 myTime = startTime
282 myIter = nIter0
283 #ifdef ALLOW_ADMTLM
284 DO k=1,4*Nr+1
285 DO j=1,sNy
286 DO i=1,sNx
287 g_objf_state_final(i,j,1,1,k) = 0.
288 ENDDO
289 ENDDO
290 ENDDO
291 #else
292 g_fc = 0.
293 #endif
294
295 CALL G_THE_MAIN_LOOP( myTime, myIter, myThid )
296 #ifdef ALLOW_ADMTLM
297 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
298 #else
299 ftlxxmemo = g_fc
300 #endif
301
302 C-- 3. reset control vector
303 CALL GRDCHK_SETXX( icvrec, TANGENT_SIMULATION,
304 & itile, jtile, layer,
305 & itilepos, jtilepos,
306 & xxmemo_ref, ierr, myThid )
307
308 #endif /* ALLOW_TANGENTLINEAR_RUN */
309
310 C******************************************************
311 C-- (C): Get gradient via finite difference perturbation
312 C******************************************************
313
314 C-- (C-1) positive perturbation:
315 localEps = ABS(grdchk_eps)
316
317 C-- get control vector component from file
318 C-- perturb it and write back to file
319 CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
320 & itile, jtile, layer,
321 & itilepos, jtilepos,
322 & xxmemo_ref, xxmemo_pert, localEps,
323 & ierr, myThid )
324
325 C-- forward run with perturbed control vector
326 myTime = startTime
327 myIter = nIter0
328 #ifdef ALLOW_AUTODIFF_OPENAD
329 CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
330 #else
331 CALL THE_MAIN_LOOP( myTime, myIter, myThid )
332 #endif
333
334 #if (defined (ALLOW_ADMTLM))
335 fcpertplus = objf_state_final(idep,jdep,1,1,1)
336 #elif (defined (ALLOW_AUTODIFF_OPENAD))
337 fcpertplus = fc%v
338 #else
339 fcpertplus = fc
340 #endif
341 WRITE(msgBuf,'(A,1PE22.14)')
342 & 'grdchk perturb(+)fc: fcpertplus =', fcpertplus
343 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
344
345 C-- Reset control vector.
346 CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
347 & itile, jtile, layer,
348 & itilepos, jtilepos,
349 & xxmemo_ref, ierr, myThid )
350
351 fcpertminus = fcref
352 IF ( useCentralDiff ) THEN
353 C-- (C-2) repeat the proceedure for a negative perturbation:
354 localEps = - ABS(grdchk_eps)
355
356 C-- get control vector component from file
357 C-- perturb it and write back to file
358 CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
359 & itile, jtile, layer,
360 & itilepos, jtilepos,
361 & xxmemo_ref, xxmemo_pert, localEps,
362 & ierr, myThid )
363
364 C-- forward run with perturbed control vector
365 myTime = startTime
366 myIter = nIter0
367 #ifdef ALLOW_AUTODIFF_OPENAD
368 CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
369 #else
370 CALL THE_MAIN_LOOP( myTime, myIter, myThid )
371 #endif
372
373 #if (defined (ALLOW_ADMTLM))
374 fcpertminus = objf_state_final(idep,jdep,1,1,1)
375 #elif (defined (ALLOW_AUTODIFF_OPENAD))
376 fcpertminus = fc%v
377 #else
378 fcpertminus = fc
379 #endif
380 WRITE(msgBuf,'(A,1PE22.14)')
381 & 'grdchk perturb(-)fc: fcpertminus =', fcpertminus
382 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
383
384 C-- Reset control vector.
385 CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
386 & itile, jtile, layer,
387 & itilepos, jtilepos,
388 & xxmemo_ref, ierr, myThid )
389
390 C-- end of if useCentralDiff ...
391 ENDIF
392
393 C******************************************************
394 C-- (D): calculate relative differences between gradients
395 C******************************************************
396
397 IF ( grdchk_eps .EQ. 0. ) THEN
398 gfd = (fcpertplus-fcpertminus)
399 ELSE
400 gfd = (fcpertplus-fcpertminus) /(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. ierr .EQ. 0 ) THEN
416 tmpplot1(itilepos,jtilepos,layer,itile,jtile) = gfd
417 tmpplot2(itilepos,jtilepos,layer,itile,jtile) = ratio_ad
418 tmpplot3(itilepos,jtilepos,layer,itile,jtile) = ratio_ftl
419 ENDIF
420
421 IF ( ierr .EQ. 0 ) THEN
422 fcrmem ( ichknum ) = fcref
423 fcppmem ( ichknum ) = fcpertplus
424 fcpmmem ( ichknum ) = fcpertminus
425 xxmemref ( ichknum ) = xxmemo_ref
426 xxmempert ( ichknum ) = xxmemo_pert
427 gfdmem ( ichknum ) = gfd
428 adxxmem ( ichknum ) = adxxmemo
429 ftlxxmem ( ichknum ) = ftlxxmemo
430 ratioadmem ( ichknum ) = ratio_ad
431 ratioftlmem ( ichknum ) = ratio_ftl
432
433 irecmem ( ichknum ) = icvrec
434 bimem ( ichknum ) = itile
435 bjmem ( ichknum ) = jtile
436 ilocmem ( ichknum ) = itilepos
437 jlocmem ( ichknum ) = jtilepos
438 klocmem ( ichknum ) = layer
439 iobcsmem ( ichknum ) = obcspos
440 icompmem ( ichknum ) = icomp
441 ichkmem ( ichknum ) = ichknum
442 itestmem ( ichknum ) = itest
443 ierrmem ( ichknum ) = ierr
444 icglomem ( ichknum ) = icglo
445 ENDIF
446
447 IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
448 WRITE(standardmessageunit,'(A)')
449 & 'grad-res -------------------------------'
450 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
451 & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
452 & layer,itile,jtile,obcspos,
453 & fcref, fcpertplus, fcpertminus
454 #ifdef ALLOW_TANGENTLINEAR_RUN
455 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
456 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
457 & icompmem(ichknum),itestmem(ichknum),
458 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
459 & ftlxxmemo, gfd, ratio_ftl
460 #else
461 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
462 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
463 & icompmem(ichknum),itestmem(ichknum),
464 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
465 & adxxmemo, gfd, ratio_ad
466 #endif
467 ENDIF
468 #ifdef ALLOW_TANGENTLINEAR_RUN
469 WRITE(msgBuf,'(A30,1PE22.14)')
470 & ' TLM ref_cost_function =', fcref
471 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
472 WRITE(msgBuf,'(A30,1PE22.14)')
473 & ' TLM tangent-lin_grad =', ftlxxmemo
474 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
475 WRITE(msgBuf,'(A30,1PE22.14)')
476 & ' TLM finite-diff_grad =', gfd
477 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
478 #else
479 WRITE(msgBuf,'(A30,1PE22.14)')
480 & ' ADM ref_cost_function =', fcref
481 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
482 WRITE(msgBuf,'(A30,1PE22.14)')
483 & ' ADM adjoint_gradient =', adxxmemo
484 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
485 WRITE(msgBuf,'(A30,1PE22.14)')
486 & ' ADM finite-diff_grad =', gfd
487 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
488 #endif
489
490 WRITE(msgBuf,'(A,I4,A,I3,A)')
491 & '====== End of gradient-check number', ichknum,
492 & ' (ierr=', ierr, ') ======='
493 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
494
495 C-- else of if ichknum < ...
496 ELSE
497 ierr_grdchk = -1
498
499 C-- end of if ichknum < ...
500 ENDIF
501
502 C-- end of do icomp = ...
503 ENDDO
504
505 IF (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) THEN
506 CALL WRITE_REC_XYZ_RL(
507 & 'grd_findiff' , tmpplot1, 1, 0, myThid)
508 CALL WRITE_REC_XYZ_RL(
509 & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
510 CALL WRITE_REC_XYZ_RL(
511 & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
512 ENDIF
513
514 C-- Print the results of the gradient check.
515 CALL GRDCHK_PRINT( ichknum, ierr_grdchk, myThid )
516
517 #endif /* ALLOW_GRDCHK */
518
519 RETURN
520 END

  ViewVC Help
Powered by ViewVC 1.1.22