/[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.36 - (show annotations) (download)
Fri Jul 6 23:10:28 2012 UTC (11 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63q, checkpoint63r
Changes since 1.35: +2 -1 lines
add explicitly #include "AD_CONFIG.h" in each fortran src file where it is needed

1 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.35 2012/02/09 17:55:54 jmc 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 #ifdef ALLOW_ADMTLM
163 fcref = objf_state_final(idep,jdep,1,1,1)
164 #else
165 fcref = fc
166 #endif
167
168 print *, 'ph-check fcref = ', fcref
169
170 do bj = jtlo, jthi
171 do bi = itlo, ithi
172 do k = 1, nr
173 do j = jmin, jmax
174 do i = imin, imax
175 tmpplot1(i,j,k,bi,bj) = 0. _d 0
176 tmpplot2(i,j,k,bi,bj) = 0. _d 0
177 tmpplot3(i,j,k,bi,bj) = 0. _d 0
178 end do
179 end do
180 end do
181 end do
182 end do
183
184 if ( useCentralDiff ) then
185 grdchk_epsfac = 2. _d 0
186 else
187 grdchk_epsfac = 1. _d 0
188 end if
189
190 WRITE(standardmessageunit,'(A)')
191 & 'grad-res -------------------------------'
192 WRITE(standardmessageunit,'(2a)')
193 & ' grad-res proc # i j k bi bj iobc',
194 & ' fc ref fc + eps fc - eps'
195 #ifdef ALLOW_TANGENTLINEAR_RUN
196 WRITE(standardmessageunit,'(2a)')
197 & ' grad-res proc # i j k bi bj iobc',
198 & ' tlm grad fd grad 1 - fd/tlm'
199 #else
200 WRITE(standardmessageunit,'(2a)')
201 & ' grad-res proc # i j k bi bj iobc',
202 & ' adj grad fd grad 1 - fd/adj'
203 #endif
204
205 c-- Compute the finite difference approximations.
206 c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
207 c-- gradient checks.
208
209 if ( nbeg .EQ. 0 )
210 & call grdchk_get_position( mythid )
211
212 do icomp = nbeg, nend, nstep
213
214 ichknum = (icomp - nbeg)/nstep + 1
215 adxxmemo = 0.
216
217 cph(
218 cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
219 cph-print & nbeg, icomp, ichknum
220 cph)
221 if (ichknum .le. maxgrdchecks ) then
222
223 c-- Determine the location of icomp on the grid.
224 if ( myProcId .EQ. grdchkwhichproc ) then
225 call grdchk_loc( icomp, ichknum,
226 & icvrec, itile, jtile, layer, obcspos,
227 & itilepos, jtilepos, icglo, itest, ierr,
228 & mythid )
229 cph(
230 cph-print print *, 'ph-grd ----- back from loc -----',
231 cph-print & icvrec, itilepos, jtilepos, layer, obcspos
232 cph)
233 else
234 icvrec = 0
235 itile = 0
236 jtile = 0
237 layer = 0
238 obcspos = 0
239 itilepos = 0
240 jtilepos = 0
241 icglo = 0
242 itest = 0
243 endif
244 _BARRIER
245
246 c******************************************************
247 c-- (A): get gradient component calculated via adjoint
248 c******************************************************
249
250 c-- get gradient component calculated via adjoint
251 if ( myProcId .EQ. grdchkwhichproc .AND.
252 & ierr .EQ. 0 ) then
253 call grdchk_getadxx( icvrec,
254 & itile, jtile, layer,
255 & itilepos, jtilepos,
256 & adxxmemo, mythid )
257 endif
258 C-- Add a global-sum call so that all proc will get the adjoint gradient
259 _GLOBAL_SUM_RL( adxxmemo, myThid )
260 c _BARRIER
261
262 #ifdef ALLOW_TANGENTLINEAR_RUN
263 c******************************************************
264 c-- (B): Get gradient component g_fc from tangent linear run:
265 c******************************************************
266 c--
267 c-- 1. perturb control vector component: xx(i)=1.
268
269 if ( myProcId .EQ. grdchkwhichproc .AND.
270 & ierr .EQ. 0 ) then
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 & mythid )
277 else
278 xxmemo_ref = 0.
279 xxmemo_pert = 0.
280 endif
281 _BARRIER
282
283 c--
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 _BARRIER
301 #ifdef ALLOW_ADMTLM
302 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
303 #else
304 ftlxxmemo = g_fc
305 #endif
306
307 c--
308 c-- 3. reset control vector
309 if ( myProcId .EQ. grdchkwhichproc .AND.
310 & ierr .EQ. 0 ) then
311 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
312 & itile, jtile, layer,
313 & itilepos, jtilepos,
314 & xxmemo_ref, mythid )
315 endif
316 _BARRIER
317
318 #endif /* ALLOW_TANGENTLINEAR_RUN */
319
320
321 c******************************************************
322 c-- (C): Get gradient via finite difference perturbation
323 c******************************************************
324
325 c-- get control vector component from file
326 c-- perturb it and write back to file
327 c-- positive perturbation
328 localEps = abs(grdchk_eps)
329 if ( myProcId .EQ. grdchkwhichproc .AND.
330 & ierr .EQ. 0 ) then
331 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
332 & itile, jtile, layer,
333 & itilepos, jtilepos,
334 & xxmemo_ref, xxmemo_pert, localEps,
335 & mythid )
336 else
337 xxmemo_ref = 0.
338 xxmemo_pert = 0.
339 endif
340 _BARRIER
341
342 c-- forward run with perturbed control vector
343 mytime = starttime
344 myiter = niter0
345 call the_main_loop( mytime, myiter, mythid )
346 #ifdef ALLOW_ADMTLM
347 fcpertplus = objf_state_final(idep,jdep,1,1,1)
348 #else
349 fcpertplus = fc
350 #endif
351 print *, 'ph-check fcpertplus = ', fcpertplus
352 _BARRIER
353
354 c-- Reset control vector.
355 if ( myProcId .EQ. grdchkwhichproc .AND.
356 & ierr .EQ. 0 ) then
357 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
358 & itile, jtile, layer,
359 & itilepos, jtilepos,
360 & xxmemo_ref, mythid )
361 endif
362 _BARRIER
363
364 fcpertminus = fcref
365 print *, 'ph-check fcpertminus = ', fcpertminus
366
367 if ( useCentralDiff ) then
368
369 c-- get control vector component from file
370 c-- perturb it and write back to file
371 c-- repeat the proceedure for a negative perturbation
372 if ( myProcId .EQ. grdchkwhichproc .AND.
373 & ierr .EQ. 0 ) then
374 localEps = - abs(grdchk_eps)
375 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
376 & itile, jtile, layer,
377 & itilepos, jtilepos,
378 & xxmemo_ref, xxmemo_pert, localEps,
379 & mythid )
380 else
381 xxmemo_ref = 0.
382 xxmemo_pert = 0.
383 endif
384 _BARRIER
385
386 c-- forward run with perturbed control vector
387 mytime = starttime
388 myiter = niter0
389 call the_main_loop( mytime, myiter, mythid )
390 _BARRIER
391 #ifdef ALLOW_ADMTLM
392 fcpertminus = objf_state_final(idep,jdep,1,1,1)
393 #else
394 fcpertminus = fc
395 #endif
396
397 c-- Reset control vector.
398 if ( myProcId .EQ. grdchkwhichproc .AND.
399 & ierr .EQ. 0 ) then
400 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
401 & itile, jtile, layer,
402 & itilepos, jtilepos,
403 & xxmemo_ref, mythid )
404 endif
405 _BARRIER
406
407 c-- end of if useCentralDiff ...
408 end if
409
410 c******************************************************
411 c-- (D): calculate relative differences between gradients
412 c******************************************************
413
414 if ( grdchk_eps .eq. 0. ) then
415 gfd = (fcpertplus-fcpertminus)
416 else
417 gfd = (fcpertplus-fcpertminus)
418 & /(grdchk_epsfac*grdchk_eps)
419 endif
420
421 if ( adxxmemo .eq. 0. ) then
422 ratio_ad = abs( adxxmemo - gfd )
423 else
424 ratio_ad = 1. - gfd/adxxmemo
425 endif
426
427 if ( ftlxxmemo .eq. 0. ) then
428 ratio_ftl = abs( ftlxxmemo - gfd )
429 else
430 ratio_ftl = 1. - gfd/ftlxxmemo
431 endif
432
433 if ( myProcId .EQ. grdchkwhichproc .AND.
434 & ierr .EQ. 0 ) then
435 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
436 & = gfd
437 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
438 & = ratio_ad
439 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
440 & = ratio_ftl
441 endif
442
443 if ( ierr .EQ. 0 ) then
444 fcrmem ( ichknum ) = fcref
445 fcppmem ( ichknum ) = fcpertplus
446 fcpmmem ( ichknum ) = fcpertminus
447 xxmemref ( ichknum ) = xxmemo_ref
448 xxmempert ( ichknum ) = xxmemo_pert
449 gfdmem ( ichknum ) = gfd
450 adxxmem ( ichknum ) = adxxmemo
451 ftlxxmem ( ichknum ) = ftlxxmemo
452 ratioadmem ( ichknum ) = ratio_ad
453 ratioftlmem ( ichknum ) = ratio_ftl
454
455 irecmem ( ichknum ) = icvrec
456 bimem ( ichknum ) = itile
457 bjmem ( ichknum ) = jtile
458 ilocmem ( ichknum ) = itilepos
459 jlocmem ( ichknum ) = jtilepos
460 klocmem ( ichknum ) = layer
461 iobcsmem ( ichknum ) = obcspos
462 icompmem ( ichknum ) = icomp
463 ichkmem ( ichknum ) = ichknum
464 itestmem ( ichknum ) = itest
465 ierrmem ( ichknum ) = ierr
466 icglomem ( ichknum ) = icglo
467 endif
468
469 if ( myProcId .EQ. grdchkwhichproc .AND.
470 & ierr .EQ. 0 ) then
471
472 WRITE(standardmessageunit,'(A)')
473 & 'grad-res -------------------------------'
474 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
475 & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
476 & layer,itile,jtile,obcspos,
477 & fcref, fcpertplus, fcpertminus
478 #ifdef ALLOW_TANGENTLINEAR_RUN
479 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
480 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
481 & icompmem(ichknum),itestmem(ichknum),
482 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
483 & ftlxxmemo, gfd, ratio_ftl
484 #else
485 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
486 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
487 & icompmem(ichknum),itestmem(ichknum),
488 & bimem(ichknum),bjmem(ichknum),obcspos,
489 & adxxmemo, gfd, ratio_ad
490 #endif
491 endif
492 #ifdef ALLOW_TANGENTLINEAR_RUN
493 WRITE(msgBuf,'(A34,1PE24.14)')
494 & ' TLM precision_derivative_cost =', fcref
495 CALL PRINT_MESSAGE
496 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
497 WRITE(msgBuf,'(A34,1PE24.14)')
498 & ' TLM precision_derivative_grad =', ftlxxmemo
499 CALL PRINT_MESSAGE
500 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
501 #else
502 WRITE(msgBuf,'(A30,1PE22.14)')
503 & ' ADM ref_cost_function =', fcref
504 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
505 & SQUEEZE_RIGHT, myThid )
506 WRITE(msgBuf,'(A30,1PE22.14)')
507 & ' ADM adjoint_gradient =', adxxmemo
508 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
509 & SQUEEZE_RIGHT, myThid )
510 WRITE(msgBuf,'(A30,1PE22.14)')
511 & ' ADM finite-diff_grad =', gfd
512 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
513 & SQUEEZE_RIGHT, myThid )
514 #endif
515
516 print *, 'ph-grd ierr ---------------------------'
517 print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
518 & ', ichknum = ', ichknum
519
520 _BARRIER
521
522 c-- else of if ( ichknum ...
523 else
524 ierr_grdchk = -1
525
526 c-- end of if ( ichknum ...
527 endif
528
529 c-- end of do icomp = ...
530 enddo
531
532 if ( myProcId .EQ. grdchkwhichproc ) then
533 CALL WRITE_REC_XYZ_RL(
534 & 'grd_findiff' , tmpplot1, 1, 0, myThid)
535 CALL WRITE_REC_XYZ_RL(
536 & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
537 CALL WRITE_REC_XYZ_RL(
538 & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
539 endif
540
541 c-- Everyone has to wait for the component to be reset.
542 _BARRIER
543
544 c-- Print the results of the gradient check.
545 call grdchk_print( ichknum, ierr_grdchk, mythid )
546
547 #endif /* ALLOW_GRDCHK */
548
549 return
550 end

  ViewVC Help
Powered by ViewVC 1.1.22