/[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.34 - (show annotations) (download)
Wed Feb 8 03:14:13 2012 UTC (12 years, 3 months ago) by jmc
Branch: MAIN
Changes since 1.33: +22 -2 lines
- all procs print a meaningfulgradient-check report
- print also the RMS of all grdcheck ratio

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

  ViewVC Help
Powered by ViewVC 1.1.22