/[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.39 - (show annotations) (download)
Tue Aug 21 19:49:41 2012 UTC (11 years, 8 months ago) by gforget
Branch: MAIN
Changes since 1.38: +2 -1 lines
- grdchk_main.F : added CALL GLOBAL_SUM_INT( ierr , myThid )
- grdchk_getadxx.F etc. : omit I/O if ierr.EQ.0

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

  ViewVC Help
Powered by ViewVC 1.1.22