/[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.37 - (show annotations) (download)
Wed Aug 15 23:05:48 2012 UTC (11 years, 9 months ago) by jmc
Branch: MAIN
Changes since 1.36: +15 -24 lines
- change format of TLM output (for testreport)
- remove all BARRIER calls (most were not needed; to check in multi-threads)

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

  ViewVC Help
Powered by ViewVC 1.1.22