/[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.32 - (show annotations) (download)
Fri Aug 12 18:36:49 2011 UTC (12 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63b
Changes since 1.31: +26 -30 lines
send adjoint gradient to all procs so that all procs print it to STDOUT
 (a hack to get gradient-check tested with testreport)

1 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.31 2011/05/24 22:40:30 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 endif
233 _BARRIER
234
235 c******************************************************
236 c-- (A): get gradient component calculated via adjoint
237 c******************************************************
238
239 c-- get gradient component calculated via adjoint
240 if ( myProcId .EQ. grdchkwhichproc .AND.
241 & ierr .EQ. 0 ) then
242 call grdchk_getadxx( icvrec,
243 & itile, jtile, layer,
244 & itilepos, jtilepos,
245 & adxxmemo, mythid )
246 endif
247 C-- Add a global-sum call so that all proc will get the adjoint gradient
248 _GLOBAL_SUM_RL( adxxmemo, myThid )
249 c _BARRIER
250
251 #ifdef ALLOW_TANGENTLINEAR_RUN
252 c******************************************************
253 c-- (B): Get gradient component g_fc from tangent linear run:
254 c******************************************************
255 c--
256 c-- 1. perturb control vector component: xx(i)=1.
257
258 if ( myProcId .EQ. grdchkwhichproc .AND.
259 & ierr .EQ. 0 ) then
260 localEps = 1. _d 0
261 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
262 & itile, jtile, layer,
263 & itilepos, jtilepos,
264 & xxmemo_ref, xxmemo_pert, localEps,
265 & mythid )
266 endif
267 _BARRIER
268
269 c--
270 c-- 2. perform tangent linear run
271 mytime = starttime
272 myiter = niter0
273 #ifdef ALLOW_ADMTLM
274 do k=1,4*Nr+1
275 do j=1,sny
276 do i=1,snx
277 g_objf_state_final(i,j,1,1,k) = 0.
278 enddo
279 enddo
280 enddo
281 #else
282 g_fc = 0.
283 #endif
284
285 call g_the_main_loop( mytime, myiter, mythid )
286 _BARRIER
287 #ifdef ALLOW_ADMTLM
288 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
289 #else
290 ftlxxmemo = g_fc
291 #endif
292
293 c--
294 c-- 3. reset control vector
295 if ( myProcId .EQ. grdchkwhichproc .AND.
296 & ierr .EQ. 0 ) then
297 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
298 & itile, jtile, layer,
299 & itilepos, jtilepos,
300 & xxmemo_ref, mythid )
301 endif
302 _BARRIER
303
304 #endif /* ALLOW_TANGENTLINEAR_RUN */
305
306
307 c******************************************************
308 c-- (C): Get gradient via finite difference perturbation
309 c******************************************************
310
311 c-- get control vector component from file
312 c-- perturb it and write back to file
313 c-- positive perturbation
314 localEps = abs(grdchk_eps)
315 if ( myProcId .EQ. grdchkwhichproc .AND.
316 & ierr .EQ. 0 ) then
317 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
318 & itile, jtile, layer,
319 & itilepos, jtilepos,
320 & xxmemo_ref, xxmemo_pert, localEps,
321 & mythid )
322 endif
323 _BARRIER
324
325 c-- forward run with perturbed control vector
326 mytime = starttime
327 myiter = niter0
328 call the_main_loop( mytime, myiter, mythid )
329 #ifdef ALLOW_ADMTLM
330 fcpertplus = objf_state_final(idep,jdep,1,1,1)
331 #else
332 fcpertplus = fc
333 #endif
334 print *, 'ph-check fcpertplus = ', fcpertplus
335 _BARRIER
336
337 c-- Reset control vector.
338 if ( myProcId .EQ. grdchkwhichproc .AND.
339 & ierr .EQ. 0 ) then
340 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
341 & itile, jtile, layer,
342 & itilepos, jtilepos,
343 & xxmemo_ref, mythid )
344 endif
345 _BARRIER
346
347 fcpertminus = fcref
348 print *, 'ph-check fcpertminus = ', fcpertminus
349
350 if ( useCentralDiff ) then
351
352 c-- get control vector component from file
353 c-- perturb it and write back to file
354 c-- repeat the proceedure for a negative perturbation
355 if ( myProcId .EQ. grdchkwhichproc .AND.
356 & ierr .EQ. 0 ) then
357 localEps = - abs(grdchk_eps)
358 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
359 & itile, jtile, layer,
360 & itilepos, jtilepos,
361 & xxmemo_ref, xxmemo_pert, localEps,
362 & mythid )
363 endif
364 _BARRIER
365
366 c-- forward run with perturbed control vector
367 mytime = starttime
368 myiter = niter0
369 call the_main_loop( mytime, myiter, mythid )
370 _BARRIER
371 #ifdef ALLOW_ADMTLM
372 fcpertminus = objf_state_final(idep,jdep,1,1,1)
373 #else
374 fcpertminus = fc
375 #endif
376
377 c-- Reset control vector.
378 if ( myProcId .EQ. grdchkwhichproc .AND.
379 & ierr .EQ. 0 ) then
380 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
381 & itile, jtile, layer,
382 & itilepos, jtilepos,
383 & xxmemo_ref, mythid )
384 endif
385 _BARRIER
386
387 c-- end of if useCentralDiff ...
388 end if
389
390 c******************************************************
391 c-- (D): calculate relative differences between gradients
392 c******************************************************
393
394 if ( myProcId .EQ. grdchkwhichproc .AND.
395 & ierr .EQ. 0 ) then
396
397 if ( grdchk_eps .eq. 0. ) then
398 gfd = (fcpertplus-fcpertminus)
399 else
400 gfd = (fcpertplus-fcpertminus)
401 & /(grdchk_epsfac*grdchk_eps)
402 endif
403
404 if ( adxxmemo .eq. 0. ) then
405 ratio_ad = abs( adxxmemo - gfd )
406 else
407 ratio_ad = 1. - gfd/adxxmemo
408 endif
409
410 if ( ftlxxmemo .eq. 0. ) then
411 ratio_ftl = abs( ftlxxmemo - gfd )
412 else
413 ratio_ftl = 1. - gfd/ftlxxmemo
414 endif
415
416 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
417 & = gfd
418 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
419 & = ratio_ad
420 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
421 & = ratio_ftl
422
423 fcrmem ( ichknum ) = fcref
424 fcppmem ( ichknum ) = fcpertplus
425 fcpmmem ( ichknum ) = fcpertminus
426 xxmemref ( ichknum ) = xxmemo_ref
427 xxmempert ( ichknum ) = xxmemo_pert
428 gfdmem ( ichknum ) = gfd
429 adxxmem ( ichknum ) = adxxmemo
430 ftlxxmem ( ichknum ) = ftlxxmemo
431 ratioadmem ( ichknum ) = ratio_ad
432 ratioftlmem ( ichknum ) = ratio_ftl
433
434 irecmem ( ichknum ) = icvrec
435 bimem ( ichknum ) = itile
436 bjmem ( ichknum ) = jtile
437 ilocmem ( ichknum ) = itilepos
438 jlocmem ( ichknum ) = jtilepos
439 klocmem ( ichknum ) = layer
440 iobcsmem ( ichknum ) = obcspos
441 icompmem ( ichknum ) = icomp
442 ichkmem ( ichknum ) = ichknum
443 itestmem ( ichknum ) = itest
444 ierrmem ( ichknum ) = ierr
445 icglomem ( ichknum ) = icglo
446
447 WRITE(standardmessageunit,'(A)')
448 & 'grad-res -------------------------------'
449 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
450 & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
451 & layer,itile,jtile,obcspos,
452 & fcref, fcpertplus, fcpertminus
453 #ifdef ALLOW_TANGENTLINEAR_RUN
454 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
455 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
456 & icompmem(ichknum),itestmem(ichknum),
457 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
458 & ftlxxmemo, gfd, ratio_ftl
459 #else
460 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
461 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
462 & icompmem(ichknum),itestmem(ichknum),
463 & bimem(ichknum),bjmem(ichknum),obcspos,
464 & adxxmemo, gfd, ratio_ad
465 #endif
466 endif
467 #ifdef ALLOW_TANGENTLINEAR_RUN
468 WRITE(msgBuf,'(A34,1PE24.14)')
469 & ' TLM precision_derivative_cost =', fcref
470 CALL PRINT_MESSAGE
471 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
472 WRITE(msgBuf,'(A34,1PE24.14)')
473 & ' TLM precision_derivative_grad =', ftlxxmemo
474 CALL PRINT_MESSAGE
475 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
476 #else
477 WRITE(msgBuf,'(A34,1PE24.14)')
478 & ' ADM precision_derivative_cost =', fcref
479 CALL PRINT_MESSAGE
480 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
481 WRITE(msgBuf,'(A34,1PE24.14)')
482 & ' ADM precision_derivative_grad =', adxxmemo
483 CALL PRINT_MESSAGE
484 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
485 #endif
486
487 print *, 'ph-grd ierr ---------------------------'
488 print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
489 & ', ichknum = ', ichknum
490
491 _BARRIER
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 ) 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-- Everyone has to wait for the component to be reset.
513 _BARRIER
514
515 c-- Print the results of the gradient check.
516 call grdchk_print( ichknum, ierr_grdchk, mythid )
517
518 #endif /* ALLOW_GRDCHK */
519
520 return
521 end

  ViewVC Help
Powered by ViewVC 1.1.22