/[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.43 - (show annotations) (download)
Fri Apr 4 21:39:04 2014 UTC (10 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64w, checkpoint64v, HEAD
Changes since 1.42: +12 -6 lines
- Start to include explicitly AUTODIFF_OPTIONS.h, COST_OPTIONS.h,
  and CTRL_OPTIONS.h in src files (to enable to skip the ECCO_CPPOPTIONS.h)
  For now, only in pkgs used in verification/hs94.1x64x5.
- Replace ALLOW_AUTODIFF_OPENAD by ALLOW_OPENAD:
  since ALLOW_OPENAD is defined in PACKAGES_CONFIG.h (any time pkg/openad
  is compiled), this simplifies/reduces which *_OPTIONS.h file to include.

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

  ViewVC Help
Powered by ViewVC 1.1.22