/[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.24 - (show annotations) (download)
Fri Sep 14 20:40:20 2007 UTC (16 years, 8 months ago) by heimbach
Branch: MAIN
Changes since 1.23: +19 -1 lines
Adding new tags in separate linesfor testreport to change AD testing:
precision_derivative_cost
precision_derivative_grad

1 C
2 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.23 2007/08/04 18:05:23 heimbach Exp $
3 C $Name: $
4
5 #include "AD_CONFIG.h"
6 #include "CPP_OPTIONS.h"
7
8 CBOI
9 C
10 C !TITLE: GRADIENT CHECK
11 C !AUTHORS: mitgcm developers ( support@mitgcm.org )
12 C !AFFILIATION: Massachussetts Institute of Technology
13 C !DATE:
14 C !INTRODUCTION: gradient check package
15 c \bv
16 c Compare the gradients calculated by the adjoint model with
17 c finite difference approximations.
18 c
19 C !CALLING SEQUENCE:
20 c
21 c the_model_main
22 c |
23 c |-- ctrl_unpack
24 c |-- adthe_main_loop - unperturbed cost function and
25 c |-- ctrl_pack adjoint gradient are computed here
26 c |
27 c |-- grdchk_main
28 c |
29 c |-- grdchk_init
30 c |-- do icomp=... - loop over control vector elements
31 c |
32 c |-- grdchk_loc - determine location of icomp on grid
33 c |
34 c |-- grdchk_getxx - get control vector component from file
35 c | perturb it and write back to file
36 c |-- grdchk_getadxx - get gradient component calculated
37 c | via adjoint
38 c |-- the_main_loop - forward run and cost evaluation
39 c | with perturbed control vector element
40 c |-- calculate ratio of adj. vs. finite difference gradient
41 c |
42 c |-- grdchk_setxx - Reset control vector element
43 c |
44 c |-- grdchk_print - print results
45 c \ev
46 CEOI
47
48 CBOP
49 C !ROUTINE: grdchk_main
50 C !INTERFACE:
51 subroutine grdchk_main( mythid )
52
53 C !DESCRIPTION: \bv
54 c ==================================================================
55 c SUBROUTINE grdchk_main
56 c ==================================================================
57 c o Compare the gradients calculated by the adjoint model with
58 c finite difference approximations.
59 c started: Christian Eckert eckert@mit.edu 24-Feb-2000
60 c continued&finished: heimbach@mit.edu: 13-Jun-2001
61 c changed: mlosch@ocean.mit.edu: 09-May-2002
62 c - added centered difference vs. 1-sided difference option
63 c - improved output format for readability
64 c - added control variable hFacC
65 c heimbach@mit.edu 24-Feb-2003
66 c - added tangent linear gradient checks
67 c - fixes for multiproc. gradient checks
68 c - added more control variables
69 c
70 c ==================================================================
71 c SUBROUTINE grdchk_main
72 c ==================================================================
73 C \ev
74
75 C !USES:
76 implicit none
77
78 c == global variables ==
79 #include "SIZE.h"
80 #include "EEPARAMS.h"
81 #include "PARAMS.h"
82 #include "grdchk.h"
83 #include "cost.h"
84 #include "ctrl.h"
85 #ifdef ALLOW_TANGENTLINEAR_RUN
86 #include "g_cost.h"
87 #endif
88
89 C !INPUT/OUTPUT PARAMETERS:
90 c == routine arguments ==
91 integer mythid
92
93 #ifdef ALLOW_GRDCHK
94 C !LOCAL VARIABLES:
95 c == local variables ==
96 integer myiter
97 _RL mytime
98 integer bi, itlo, ithi
99 integer bj, jtlo, jthi
100 integer i, imin, imax
101 integer j, jmin, jmax
102 integer k
103
104 integer icomp
105 integer ichknum
106 integer icvrec
107 integer jtile
108 integer itile
109 integer layer
110 integer obcspos
111 integer itilepos
112 integer jtilepos
113 integer icglo
114 integer itest
115 integer ierr
116 integer ierr_grdchk
117 _RL gfd
118 _RL fcref
119 _RL fcpertplus, fcpertminus
120 _RL ratio_ad
121 _RL ratio_ftl
122 _RL xxmemo_ref
123 _RL xxmemo_pert
124 _RL adxxmemo
125 _RL ftlxxmemo
126 _RL localEps
127 _RL grdchk_epsfac
128
129 _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130 _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
131 _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
132
133 CHARACTER*(MAX_LEN_MBUF) msgBuf
134
135 c == end of interface ==
136 CEOP
137
138 c-- Set the loop ranges.
139 jtlo = mybylo(mythid)
140 jthi = mybyhi(mythid)
141 itlo = mybxlo(mythid)
142 ithi = mybxhi(mythid)
143 jmin = 1
144 jmax = sny
145 imin = 1
146 imax = snx
147
148 print *, 'ph-check entering grdchk_main '
149
150 c-- initialise variables
151 call grdchk_init( mythid )
152
153 c-- Compute the adjoint models' gradients.
154 c-- Compute the unperturbed cost function.
155 cph Gradient via adjoint has already been computed,
156 cph and so has unperturbed cost function,
157 cph assuming all xx_ fields are initialised to zero.
158
159 ierr_grdchk = 0
160 #ifdef ALLOW_ADMTLM
161 fcref = objf_state_final(idep,jdep,1,1,1)
162 #else
163 fcref = fc
164 #endif
165
166 print *, 'ph-check fcref = ', fcref
167
168 do bj = jtlo, jthi
169 do bi = itlo, ithi
170 do k = 1, nr
171 do j = jmin, jmax
172 do i = imin, imax
173 tmpplot1(i,j,k,bi,bj) = 0. _d 0
174 tmpplot2(i,j,k,bi,bj) = 0. _d 0
175 tmpplot3(i,j,k,bi,bj) = 0. _d 0
176 end do
177 end do
178 end do
179 end do
180 end do
181
182 if ( useCentralDiff ) then
183 grdchk_epsfac = 2. _d 0
184 else
185 grdchk_epsfac = 1. _d 0
186 end if
187
188 WRITE(standardmessageunit,'(A)')
189 & 'grad-res -------------------------------'
190 WRITE(standardmessageunit,'(2a)')
191 & ' grad-res proc # i j k bi bj iobc',
192 & ' fc ref fc + eps fc - eps'
193 #ifdef ALLOW_TANGENTLINEAR_RUN
194 WRITE(standardmessageunit,'(2a)')
195 & ' grad-res proc # i j k bi bj iobc',
196 & ' tlm grad fd grad 1 - fd/tlm'
197 #else
198 WRITE(standardmessageunit,'(2a)')
199 & ' grad-res proc # i j k bi bj iobc',
200 & ' adj grad fd grad 1 - fd/adj'
201 #endif
202
203 c-- Compute the finite difference approximations.
204 c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
205 c-- gradient checks.
206
207 if ( nbeg .EQ. 0 )
208 & call grdchk_get_position( mythid )
209
210 do icomp = nbeg, nend, nstep
211
212 ichknum = (icomp - nbeg)/nstep + 1
213
214 cph(
215 cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
216 cph-print & nbeg, icomp, ichknum
217 cph)
218 if (ichknum .le. maxgrdchecks ) then
219
220 c-- Determine the location of icomp on the grid.
221 if ( myProcId .EQ. grdchkwhichproc ) then
222 call grdchk_loc( icomp, ichknum,
223 & icvrec, itile, jtile, layer, obcspos,
224 & itilepos, jtilepos, icglo, itest, ierr,
225 & mythid )
226 cph(
227 cph-print print *, 'ph-grd ----- back from loc -----',
228 cph-print & icvrec, itilepos, jtilepos, layer, obcspos
229 cph)
230 endif
231 _BARRIER
232
233 c******************************************************
234 c-- (A): get gradient component calculated via adjoint
235 c******************************************************
236
237 c-- get gradient component calculated via adjoint
238 if ( myProcId .EQ. grdchkwhichproc .AND.
239 & ierr .EQ. 0 ) then
240 call grdchk_getadxx( icvrec,
241 & itile, jtile, layer,
242 & itilepos, jtilepos,
243 & adxxmemo, mythid )
244 endif
245 _BARRIER
246
247 #ifdef ALLOW_TANGENTLINEAR_RUN
248 c******************************************************
249 c-- (B): Get gradient component g_fc from tangent linear run:
250 c******************************************************
251 c--
252 c-- 1. perturb control vector component: xx(i)=1.
253 ftlxxmemo = 0.
254
255 if ( myProcId .EQ. grdchkwhichproc .AND.
256 & ierr .EQ. 0 ) then
257 localEps = 1. _d 0
258 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
259 & itile, jtile, layer,
260 & itilepos, jtilepos,
261 & xxmemo_ref, xxmemo_pert, localEps,
262 & mythid )
263 endif
264 _BARRIER
265
266 c--
267 c-- 2. perform tangent linear run
268 mytime = starttime
269 myiter = niter0
270 #ifdef ALLOW_ADMTLM
271 do k=1,4*Nr+1
272 do j=1,sny
273 do i=1,snx
274 g_objf_state_final(i,j,1,1,k) = 0.
275 enddo
276 enddo
277 enddo
278 #else
279 g_fc = 0.
280 #endif
281
282 call g_the_main_loop( mytime, myiter, mythid )
283 _BARRIER
284 #ifdef ALLOW_ADMTLM
285 ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
286 #else
287 ftlxxmemo = g_fc
288 #endif
289
290 c--
291 c-- 3. reset control vector
292 if ( myProcId .EQ. grdchkwhichproc .AND.
293 & ierr .EQ. 0 ) then
294 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
295 & itile, jtile, layer,
296 & itilepos, jtilepos,
297 & xxmemo_ref, mythid )
298 endif
299 _BARRIER
300
301 #endif /* ALLOW_TANGENTLINEAR_RUN */
302
303
304 c******************************************************
305 c-- (C): Get gradient via finite difference perturbation
306 c******************************************************
307
308 c-- get control vector component from file
309 c-- perturb it and write back to file
310 c-- positive perturbation
311 localEps = abs(grdchk_eps)
312 if ( myProcId .EQ. grdchkwhichproc .AND.
313 & ierr .EQ. 0 ) then
314 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
315 & itile, jtile, layer,
316 & itilepos, jtilepos,
317 & xxmemo_ref, xxmemo_pert, localEps,
318 & mythid )
319 endif
320 _BARRIER
321
322 c-- forward run with perturbed control vector
323 mytime = starttime
324 myiter = niter0
325 call the_main_loop( mytime, myiter, mythid )
326 #ifdef ALLOW_ADMTLM
327 fcpertplus = objf_state_final(idep,jdep,1,1,1)
328 #else
329 fcpertplus = fc
330 #endif
331 print *, 'ph-check fcpertplus = ', fcpertplus
332 _BARRIER
333
334 c-- Reset control vector.
335 if ( myProcId .EQ. grdchkwhichproc .AND.
336 & ierr .EQ. 0 ) then
337 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
338 & itile, jtile, layer,
339 & itilepos, jtilepos,
340 & xxmemo_ref, mythid )
341 endif
342 _BARRIER
343
344 fcpertminus = fcref
345 print *, 'ph-check fcpertminus = ', fcpertminus
346
347 if ( useCentralDiff ) then
348
349 c-- get control vector component from file
350 c-- perturb it and write back to file
351 c-- repeat the proceedure for a negative perturbation
352 if ( myProcId .EQ. grdchkwhichproc .AND.
353 & ierr .EQ. 0 ) then
354 localEps = - abs(grdchk_eps)
355 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
356 & itile, jtile, layer,
357 & itilepos, jtilepos,
358 & xxmemo_ref, xxmemo_pert, localEps,
359 & mythid )
360 endif
361 _BARRIER
362
363 c-- forward run with perturbed control vector
364 mytime = starttime
365 myiter = niter0
366 call the_main_loop( mytime, myiter, mythid )
367 _BARRIER
368 #ifdef ALLOW_ADMTLM
369 fcpertminus = objf_state_final(idep,jdep,1,1,1)
370 #else
371 fcpertminus = fc
372 #endif
373
374 c-- Reset control vector.
375 if ( myProcId .EQ. grdchkwhichproc .AND.
376 & ierr .EQ. 0 ) then
377 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
378 & itile, jtile, layer,
379 & itilepos, jtilepos,
380 & xxmemo_ref, mythid )
381 endif
382 _BARRIER
383
384 c-- end of if useCentralDiff ...
385 end if
386
387 c******************************************************
388 c-- (D): calculate relative differences between gradients
389 c******************************************************
390
391 if ( myProcId .EQ. grdchkwhichproc .AND.
392 & ierr .EQ. 0 ) then
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 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
414 & = gfd
415 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
416 & = ratio_ad
417 tmpplot3(itilepos,jtilepos,layer,itile,jtile)
418 & = ratio_ftl
419
420 fcrmem ( ichknum ) = fcref
421 fcppmem ( ichknum ) = fcpertplus
422 fcpmmem ( ichknum ) = fcpertminus
423 xxmemref ( ichknum ) = xxmemo_ref
424 xxmempert ( ichknum ) = xxmemo_pert
425 gfdmem ( ichknum ) = gfd
426 adxxmem ( ichknum ) = adxxmemo
427 ftlxxmem ( ichknum ) = ftlxxmemo
428 ratioadmem ( ichknum ) = ratio_ad
429 ratioftlmem ( ichknum ) = ratio_ftl
430
431 irecmem ( ichknum ) = icvrec
432 bimem ( ichknum ) = itile
433 bjmem ( ichknum ) = jtile
434 ilocmem ( ichknum ) = itilepos
435 jlocmem ( ichknum ) = jtilepos
436 klocmem ( ichknum ) = layer
437 iobcsmem ( ichknum ) = obcspos
438 icompmem ( ichknum ) = icomp
439 ichkmem ( ichknum ) = ichknum
440 itestmem ( ichknum ) = itest
441 ierrmem ( ichknum ) = ierr
442 icglomem ( ichknum ) = icglo
443
444 WRITE(standardmessageunit,'(A)')
445 & 'grad-res -------------------------------'
446 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
447 & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
448 & layer,itile,jtile,obcspos,
449 & fcref, fcpertplus, fcpertminus
450 #ifdef ALLOW_TANGENTLINEAR_RUN
451 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
452 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
453 & icompmem(ichknum),itestmem(ichknum),
454 & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
455 & ftlxxmemo, gfd, ratio_ftl
456 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
457 & 'precision_grdchk_result TLM ', fcref, ftlxxmemo
458 CALL PRINT_MESSAGE
459 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
460 c
461 WRITE(msgBuf,'(A34,1PE24.14)')
462 & 'precision_derivative_cost TLM ', fcref
463 CALL PRINT_MESSAGE
464 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
465 WRITE(msgBuf,'(A34,1PE24.14)')
466 & 'precision_derivative_grad TLM ', ftlxxmemo
467 CALL PRINT_MESSAGE
468 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
469 #else
470 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
471 & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
472 & icompmem(ichknum),itestmem(ichknum),
473 & bimem(ichknum),bjmem(ichknum),obcspos,
474 & adxxmemo, gfd, ratio_ad
475 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
476 & 'precision_grdchk_result ADM ', fcref, adxxmemo
477 CALL PRINT_MESSAGE
478 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
479 c
480 WRITE(msgBuf,'(A34,1PE24.14)')
481 & 'precision_derivative_cost ADM ', fcref
482 CALL PRINT_MESSAGE
483 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
484 WRITE(msgBuf,'(A34,1PE24.14)')
485 & 'precision_derivative_grad ADM ', ftlxxmemo
486 CALL PRINT_MESSAGE
487 & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
488 #endif
489
490 endif
491
492 print *, 'ph-grd ierr ---------------------------'
493 print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
494 & ', ichknum = ', ichknum
495
496 _BARRIER
497
498 c-- else of if ( ichknum ...
499 else
500 ierr_grdchk = -1
501
502 c-- end of if ( ichknum ...
503 endif
504
505 c-- end of do icomp = ...
506 enddo
507
508 if ( myProcId .EQ. grdchkwhichproc ) then
509 CALL WRITE_REC_XYZ_RL(
510 & 'grd_findiff' , tmpplot1, 1, 0, myThid)
511 CALL WRITE_REC_XYZ_RL(
512 & 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
513 CALL WRITE_REC_XYZ_RL(
514 & 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
515 endif
516
517 c-- Everyone has to wait for the component to be reset.
518 _BARRIER
519
520 c-- Print the results of the gradient check.
521 call grdchk_print( ichknum, ierr_grdchk, mythid )
522
523 #endif /* ALLOW_GRDCHK */
524
525 end

  ViewVC Help
Powered by ViewVC 1.1.22