/[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.40 - (show annotations) (download)
Wed Aug 22 19:34:51 2012 UTC (11 years, 8 months ago) by heimbach
Branch: MAIN
Changes since 1.39: +20 -4 lines
Merge OpenAD specific stuff to pkg/grdchk

1 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.39 2012/08/21 19:49:41 gforget 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_getxx - get control vector component from file
34 c | perturb it and write back to file
35 c |-- grdchk_getadxx - get gradient component calculated
36 c | via adjoint
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 ==================================================================
70 c SUBROUTINE grdchk_main
71 c ==================================================================
72 C \ev
73
74 C !USES:
75 implicit none
76
77 c == global variables ==
78 #include "SIZE.h"
79 #include "EEPARAMS.h"
80 #include "PARAMS.h"
81 #include "grdchk.h"
82 #include "cost.h"
83 #include "ctrl.h"
84 #ifdef ALLOW_TANGENTLINEAR_RUN
85 #include "g_cost.h"
86 #endif
87
88 C !INPUT/OUTPUT PARAMETERS:
89 c == routine arguments ==
90 integer mythid
91
92 #ifdef ALLOW_GRDCHK
93 C !LOCAL VARIABLES:
94 c == local variables ==
95 integer myiter
96 _RL mytime
97 integer bi, itlo, ithi
98 integer bj, jtlo, jthi
99 integer i, imin, imax
100 integer j, jmin, jmax
101 integer k
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
128 _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
129 _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130 _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
131
132 CHARACTER*(MAX_LEN_MBUF) msgBuf
133
134 c == end of interface ==
135 CEOP
136
137 c-- Set the loop ranges.
138 jtlo = mybylo(mythid)
139 jthi = mybyhi(mythid)
140 itlo = mybxlo(mythid)
141 ithi = mybxhi(mythid)
142 jmin = 1
143 jmax = sny
144 imin = 1
145 imax = snx
146
147 print *, 'ph-check entering grdchk_main '
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
170 print *, 'ph-check fcref = ', fcref
171
172 do bj = jtlo, jthi
173 do bi = itlo, ithi
174 do k = 1, nr
175 do j = jmin, jmax
176 do i = imin, imax
177 tmpplot1(i,j,k,bi,bj) = 0. _d 0
178 tmpplot2(i,j,k,bi,bj) = 0. _d 0
179 tmpplot3(i,j,k,bi,bj) = 0. _d 0
180 end do
181 end do
182 end do
183 end do
184 end do
185
186 if ( useCentralDiff ) then
187 grdchk_epsfac = 2. _d 0
188 else
189 grdchk_epsfac = 1. _d 0
190 end if
191
192 WRITE(standardmessageunit,'(A)')
193 & 'grad-res -------------------------------'
194 WRITE(standardmessageunit,'(2a)')
195 & ' grad-res proc # i j k bi bj iobc',
196 & ' fc ref fc + eps fc - eps'
197 #ifdef ALLOW_TANGENTLINEAR_RUN
198 WRITE(standardmessageunit,'(2a)')
199 & ' grad-res proc # i j k bi bj iobc',
200 & ' tlm grad fd grad 1 - fd/tlm'
201 #else
202 WRITE(standardmessageunit,'(2a)')
203 & ' grad-res proc # i j k bi bj iobc',
204 & ' adj grad fd grad 1 - fd/adj'
205 #endif
206
207 c-- Compute the finite difference approximations.
208 c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
209 c-- gradient checks.
210
211 if ( nbeg .EQ. 0 )
212 & call grdchk_get_position( mythid )
213
214 do icomp = nbeg, nend, nstep
215
216 ichknum = (icomp - nbeg)/nstep + 1
217 adxxmemo = 0.
218
219 cph(
220 cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
221 cph-print & nbeg, icomp, ichknum
222 cph)
223 if (ichknum .le. maxgrdchecks ) then
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 cph(
232 cph-print print *, 'ph-grd ----- back from loc -----',
233 cph-print & icvrec, itilepos, jtilepos, layer, obcspos
234 cph)
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 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 c--
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--
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 #ifdef ALLOW_ADMTLM
296 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
297 #else
298 ftlxxmemo = g_fc
299 #endif
300
301 c--
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
311 c******************************************************
312 c-- (C): Get gradient via finite difference perturbation
313 c******************************************************
314
315 c-- get control vector component from file
316 c-- perturb it and write back to file
317 c-- positive perturbation
318 localEps = abs(grdchk_eps)
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 print *, 'ph-check fcpertplus = ', fcpertplus
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 print *, 'ph-check fcpertminus = ', fcpertminus
351
352 if ( useCentralDiff ) then
353
354 c-- get control vector component from file
355 c-- perturb it and write back to file
356 c-- repeat the proceedure for a negative perturbation
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 & 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
381 c-- Reset control vector.
382 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
383 & itile, jtile, layer,
384 & itilepos, jtilepos,
385 & xxmemo_ref, ierr, mythid )
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 ( grdchk_eps .eq. 0. ) then
395 gfd = (fcpertplus-fcpertminus)
396 else
397 gfd = (fcpertplus-fcpertminus)
398 & /(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.
414 & ierr .EQ. 0 ) then
415 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
416 & = gfd
417 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
418 & = ratio_ad
419 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
420 & = ratio_ftl
421 endif
422
423 if ( ierr .EQ. 0 ) then
424 fcrmem ( ichknum ) = fcref
425 fcppmem ( ichknum ) = fcpertplus
426 fcpmmem ( ichknum ) = fcpertminus
427 xxmemref ( ichknum ) = xxmemo_ref
428 xxmempert ( ichknum ) = xxmemo_pert
429 gfdmem ( ichknum ) = gfd
430 adxxmem ( ichknum ) = adxxmemo
431 ftlxxmem ( ichknum ) = ftlxxmemo
432 ratioadmem ( ichknum ) = ratio_ad
433 ratioftlmem ( ichknum ) = ratio_ftl
434
435 irecmem ( ichknum ) = icvrec
436 bimem ( ichknum ) = itile
437 bjmem ( ichknum ) = jtile
438 ilocmem ( ichknum ) = itilepos
439 jlocmem ( ichknum ) = jtilepos
440 klocmem ( ichknum ) = layer
441 iobcsmem ( ichknum ) = obcspos
442 icompmem ( ichknum ) = icomp
443 ichkmem ( ichknum ) = ichknum
444 itestmem ( ichknum ) = itest
445 ierrmem ( ichknum ) = ierr
446 icglomem ( ichknum ) = icglo
447 endif
448
449 if ( myProcId .EQ. grdchkwhichproc .AND.
450 & ierr .EQ. 0 ) then
451
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),obcspos,
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, standardMessageUnit,
476 & SQUEEZE_RIGHT, myThid )
477 WRITE(msgBuf,'(A30,1PE22.14)')
478 & ' TLM tangent-lin_grad =', ftlxxmemo
479 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
480 & SQUEEZE_RIGHT, myThid )
481 WRITE(msgBuf,'(A30,1PE22.14)')
482 & ' TLM finite-diff_grad =', gfd
483 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
484 & SQUEEZE_RIGHT, myThid )
485 #else
486 WRITE(msgBuf,'(A30,1PE22.14)')
487 & ' ADM ref_cost_function =', fcref
488 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
489 & SQUEEZE_RIGHT, myThid )
490 WRITE(msgBuf,'(A30,1PE22.14)')
491 & ' ADM adjoint_gradient =', adxxmemo
492 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
493 & SQUEEZE_RIGHT, myThid )
494 WRITE(msgBuf,'(A30,1PE22.14)')
495 & ' ADM finite-diff_grad =', gfd
496 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
497 & SQUEEZE_RIGHT, myThid )
498 #endif
499
500 print *, 'ph-grd ierr ---------------------------'
501 print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
502 & ', ichknum = ', ichknum
503
504 c-- else of if ( ichknum ...
505 else
506 ierr_grdchk = -1
507
508 c-- end of if ( ichknum ...
509 endif
510
511 c-- end of do icomp = ...
512 enddo
513
514 if (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) then
515 CALL WRITE_REC_XYZ_RL(
516 & 'grd_findiff' , tmpplot1, 1, 0, myThid)
517 CALL WRITE_REC_XYZ_RL(
518 & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
519 CALL WRITE_REC_XYZ_RL(
520 & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
521 endif
522
523 c-- Everyone has to wait for the component to be reset.
524 c _BARRIER
525
526 c-- Print the results of the gradient check.
527 call grdchk_print( ichknum, ierr_grdchk, mythid )
528
529 #endif /* ALLOW_GRDCHK */
530
531 return
532 end

  ViewVC Help
Powered by ViewVC 1.1.22