/[MITgcm]/MITgcm/pkg/grdchk/grdchk_main.F
ViewVC logotype

Annotation of /MITgcm/pkg/grdchk/grdchk_main.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.27 - (hide annotations) (download)
Mon Dec 3 22:38:27 2007 UTC (16 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint61f, checkpoint61n, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.26: +3 -2 lines
Always initialise adxxmemo, ftlxxmemo

1 edhill 1.10 C
2 heimbach 1.27 C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.26 2007/09/21 18:53:36 jmc Exp $
3 heimbach 1.11 C $Name: $
4 heimbach 1.2
5 edhill 1.10 #include "AD_CONFIG.h"
6 heimbach 1.2 #include "CPP_OPTIONS.h"
7    
8 heimbach 1.3 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 heimbach 1.2
53 heimbach 1.3 C !DESCRIPTION: \bv
54 heimbach 1.2 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 heimbach 1.4 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 heimbach 1.7 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 heimbach 1.4 c
70 heimbach 1.2 c ==================================================================
71     c SUBROUTINE grdchk_main
72     c ==================================================================
73 heimbach 1.3 C \ev
74 heimbach 1.2
75 heimbach 1.3 C !USES:
76 heimbach 1.2 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 heimbach 1.21 #include "ctrl.h"
85 heimbach 1.9 #ifdef ALLOW_TANGENTLINEAR_RUN
86     #include "g_cost.h"
87     #endif
88 heimbach 1.2
89 heimbach 1.3 C !INPUT/OUTPUT PARAMETERS:
90 heimbach 1.2 c == routine arguments ==
91     integer mythid
92    
93 heimbach 1.11 #ifdef ALLOW_GRDCHK
94 heimbach 1.3 C !LOCAL VARIABLES:
95 heimbach 1.2 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 heimbach 1.8 integer obcspos
111 heimbach 1.2 integer itilepos
112     integer jtilepos
113 heimbach 1.19 integer icglo
114 heimbach 1.2 integer itest
115     integer ierr
116     integer ierr_grdchk
117     _RL gfd
118     _RL fcref
119 heimbach 1.4 _RL fcpertplus, fcpertminus
120 heimbach 1.6 _RL ratio_ad
121     _RL ratio_ftl
122 heimbach 1.2 _RL xxmemo_ref
123     _RL xxmemo_pert
124     _RL adxxmemo
125 heimbach 1.6 _RL ftlxxmemo
126     _RL localEps
127     _RL grdchk_epsfac
128    
129 heimbach 1.7 _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 heimbach 1.12 CHARACTER*(MAX_LEN_MBUF) msgBuf
134    
135 heimbach 1.2 c == end of interface ==
136 heimbach 1.3 CEOP
137 heimbach 1.2
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 heimbach 1.16 print *, 'ph-check entering grdchk_main '
149    
150 heimbach 1.2 c-- initialise variables
151     call grdchk_init( mythid )
152    
153     c-- Compute the adjoint models' gradients.
154     c-- Compute the unperturbed cost function.
155 heimbach 1.7 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 heimbach 1.2
159 heimbach 1.9 ierr_grdchk = 0
160 heimbach 1.27 adxxmemo = 0.
161     ftlxxmemo = 0.
162 heimbach 1.21 #ifdef ALLOW_ADMTLM
163     fcref = objf_state_final(idep,jdep,1,1,1)
164     #else
165 heimbach 1.2 fcref = fc
166 heimbach 1.21 #endif
167 heimbach 1.9
168     print *, 'ph-check fcref = ', fcref
169 heimbach 1.2
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 heimbach 1.6 tmpplot3(i,j,k,bi,bj) = 0. _d 0
178 heimbach 1.2 end do
179     end do
180     end do
181     end do
182     end do
183    
184 heimbach 1.4 if ( useCentralDiff ) then
185     grdchk_epsfac = 2. _d 0
186     else
187     grdchk_epsfac = 1. _d 0
188     end if
189    
190 heimbach 1.23 WRITE(standardmessageunit,'(A)')
191     & 'grad-res -------------------------------'
192     WRITE(standardmessageunit,'(2a)')
193 heimbach 1.19 & ' grad-res proc # i j k bi bj iobc',
194 heimbach 1.18 & ' fc ref fc + eps fc - eps'
195 heimbach 1.7 #ifdef ALLOW_TANGENTLINEAR_RUN
196 heimbach 1.23 WRITE(standardmessageunit,'(2a)')
197 heimbach 1.19 & ' grad-res proc # i j k bi bj iobc',
198 heimbach 1.18 & ' tlm grad fd grad 1 - fd/tlm'
199 heimbach 1.7 #else
200 heimbach 1.23 WRITE(standardmessageunit,'(2a)')
201 heimbach 1.19 & ' grad-res proc # i j k bi bj iobc',
202 heimbach 1.18 & ' adj grad fd grad 1 - fd/adj'
203 heimbach 1.7 #endif
204    
205 heimbach 1.2 c-- Compute the finite difference approximations.
206     c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
207     c-- gradient checks.
208 heimbach 1.14
209 heimbach 1.17 if ( nbeg .EQ. 0 )
210     & call grdchk_get_position( mythid )
211 heimbach 1.2
212 heimbach 1.7 do icomp = nbeg, nend, nstep
213 heimbach 1.2
214 heimbach 1.7 ichknum = (icomp - nbeg)/nstep + 1
215 heimbach 1.2
216 heimbach 1.19 cph(
217     cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
218     cph-print & nbeg, icomp, ichknum
219     cph)
220 heimbach 1.7 if (ichknum .le. maxgrdchecks ) then
221 heimbach 1.2
222 heimbach 1.7 c-- Determine the location of icomp on the grid.
223     if ( myProcId .EQ. grdchkwhichproc ) then
224     call grdchk_loc( icomp, ichknum,
225 heimbach 1.8 & icvrec, itile, jtile, layer, obcspos,
226 heimbach 1.19 & itilepos, jtilepos, icglo, itest, ierr,
227 heimbach 1.7 & mythid )
228 heimbach 1.18 cph(
229 heimbach 1.19 cph-print print *, 'ph-grd ----- back from loc -----',
230     cph-print & icvrec, itilepos, jtilepos, layer, obcspos
231 heimbach 1.18 cph)
232 heimbach 1.7 endif
233     _BARRIER
234 heimbach 1.2
235 heimbach 1.6 c******************************************************
236     c-- (A): get gradient component calculated via adjoint
237     c******************************************************
238 heimbach 1.7
239     c-- get gradient component calculated via adjoint
240     if ( myProcId .EQ. grdchkwhichproc .AND.
241     & ierr .EQ. 0 ) then
242     call grdchk_getadxx( icvrec,
243     & itile, jtile, layer,
244     & itilepos, jtilepos,
245     & adxxmemo, mythid )
246     endif
247     _BARRIER
248 heimbach 1.4
249 heimbach 1.6 #ifdef ALLOW_TANGENTLINEAR_RUN
250     c******************************************************
251     c-- (B): Get gradient component g_fc from tangent linear run:
252     c******************************************************
253     c--
254     c-- 1. perturb control vector component: xx(i)=1.
255    
256 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
257     & ierr .EQ. 0 ) then
258     localEps = 1. _d 0
259     call grdchk_getxx( icvrec, TANGENT_SIMULATION,
260     & itile, jtile, layer,
261     & itilepos, jtilepos,
262     & xxmemo_ref, xxmemo_pert, localEps,
263     & mythid )
264     endif
265     _BARRIER
266 heimbach 1.6
267     c--
268     c-- 2. perform tangent linear run
269 heimbach 1.7 mytime = starttime
270     myiter = niter0
271 heimbach 1.21 #ifdef ALLOW_ADMTLM
272     do k=1,4*Nr+1
273     do j=1,sny
274     do i=1,snx
275     g_objf_state_final(i,j,1,1,k) = 0.
276     enddo
277     enddo
278     enddo
279     #else
280 heimbach 1.7 g_fc = 0.
281 heimbach 1.21 #endif
282    
283 heimbach 1.7 call g_the_main_loop( mytime, myiter, mythid )
284 heimbach 1.21 _BARRIER
285     #ifdef ALLOW_ADMTLM
286     ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
287     #else
288 heimbach 1.7 ftlxxmemo = g_fc
289 heimbach 1.21 #endif
290    
291 heimbach 1.6 c--
292     c-- 3. reset control vector
293 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
294     & ierr .EQ. 0 ) then
295     call grdchk_setxx( icvrec, TANGENT_SIMULATION,
296     & itile, jtile, layer,
297     & itilepos, jtilepos,
298     & xxmemo_ref, mythid )
299     endif
300     _BARRIER
301    
302     #endif /* ALLOW_TANGENTLINEAR_RUN */
303    
304 heimbach 1.6
305     c******************************************************
306     c-- (C): Get gradient via finite difference perturbation
307     c******************************************************
308    
309 heimbach 1.2 c-- get control vector component from file
310 heimbach 1.7 c-- perturb it and write back to file
311 heimbach 1.6 c-- positive perturbation
312 heimbach 1.7 localEps = abs(grdchk_eps)
313     if ( myProcId .EQ. grdchkwhichproc .AND.
314     & ierr .EQ. 0 ) then
315     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
316     & itile, jtile, layer,
317     & itilepos, jtilepos,
318     & xxmemo_ref, xxmemo_pert, localEps,
319     & mythid )
320     endif
321     _BARRIER
322 heimbach 1.2
323 heimbach 1.4 c-- forward run with perturbed control vector
324 heimbach 1.7 mytime = starttime
325     myiter = niter0
326     call the_main_loop( mytime, myiter, mythid )
327 heimbach 1.21 #ifdef ALLOW_ADMTLM
328     fcpertplus = objf_state_final(idep,jdep,1,1,1)
329     #else
330 heimbach 1.7 fcpertplus = fc
331 heimbach 1.21 #endif
332 heimbach 1.13 print *, 'ph-check fcpertplus = ', fcpertplus
333 heimbach 1.7 _BARRIER
334 heimbach 1.4
335     c-- Reset control vector.
336 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
337     & ierr .EQ. 0 ) then
338     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
339     & itile, jtile, layer,
340     & itilepos, jtilepos,
341     & xxmemo_ref, mythid )
342     endif
343     _BARRIER
344 heimbach 1.2
345 heimbach 1.7 fcpertminus = fcref
346 heimbach 1.13 print *, 'ph-check fcpertminus = ', fcpertminus
347 heimbach 1.4
348 heimbach 1.7 if ( useCentralDiff ) then
349 heimbach 1.4
350 heimbach 1.6 c-- get control vector component from file
351 heimbach 1.7 c-- perturb it and write back to file
352 heimbach 1.4 c-- repeat the proceedure for a negative perturbation
353 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
354     & ierr .EQ. 0 ) then
355     localEps = - abs(grdchk_eps)
356     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
357     & itile, jtile, layer,
358     & itilepos, jtilepos,
359     & xxmemo_ref, xxmemo_pert, localEps,
360     & mythid )
361     endif
362     _BARRIER
363 heimbach 1.4
364 heimbach 1.2 c-- forward run with perturbed control vector
365 heimbach 1.7 mytime = starttime
366     myiter = niter0
367     call the_main_loop( mytime, myiter, mythid )
368     _BARRIER
369 heimbach 1.21 #ifdef ALLOW_ADMTLM
370     fcpertminus = objf_state_final(idep,jdep,1,1,1)
371     #else
372 heimbach 1.7 fcpertminus = fc
373 heimbach 1.21 #endif
374 heimbach 1.2
375 heimbach 1.4 c-- Reset control vector.
376 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
377     & ierr .EQ. 0 ) then
378     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
379     & itile, jtile, layer,
380     & itilepos, jtilepos,
381     & xxmemo_ref, mythid )
382     endif
383     _BARRIER
384    
385     c-- end of if useCentralDiff ...
386     end if
387 heimbach 1.4
388 heimbach 1.6 c******************************************************
389     c-- (D): calculate relative differences between gradients
390     c******************************************************
391    
392 heimbach 1.7 if ( myProcId .EQ. grdchkwhichproc .AND.
393     & ierr .EQ. 0 ) then
394 heimbach 1.2
395 heimbach 1.7 if ( grdchk_eps .eq. 0. ) then
396     gfd = (fcpertplus-fcpertminus)
397     else
398     gfd = (fcpertplus-fcpertminus)
399     & /(grdchk_epsfac*grdchk_eps)
400     endif
401 heimbach 1.2
402 heimbach 1.7 if ( adxxmemo .eq. 0. ) then
403     ratio_ad = abs( adxxmemo - gfd )
404     else
405     ratio_ad = 1. - gfd/adxxmemo
406     endif
407    
408     if ( ftlxxmemo .eq. 0. ) then
409     ratio_ftl = abs( ftlxxmemo - gfd )
410 heimbach 1.2 else
411 heimbach 1.7 ratio_ftl = 1. - gfd/ftlxxmemo
412 heimbach 1.2 endif
413 heimbach 1.7
414     tmpplot1(itilepos,jtilepos,layer,itile,jtile)
415     & = gfd
416     tmpplot2(itilepos,jtilepos,layer,itile,jtile)
417     & = ratio_ad
418     tmpplot3(itilepos,jtilepos,layer,itile,jtile)
419     & = ratio_ftl
420    
421     fcrmem ( ichknum ) = fcref
422     fcppmem ( ichknum ) = fcpertplus
423     fcpmmem ( ichknum ) = fcpertminus
424     xxmemref ( ichknum ) = xxmemo_ref
425     xxmempert ( ichknum ) = xxmemo_pert
426     gfdmem ( ichknum ) = gfd
427     adxxmem ( ichknum ) = adxxmemo
428     ftlxxmem ( ichknum ) = ftlxxmemo
429     ratioadmem ( ichknum ) = ratio_ad
430     ratioftlmem ( ichknum ) = ratio_ftl
431    
432     irecmem ( ichknum ) = icvrec
433     bimem ( ichknum ) = itile
434     bjmem ( ichknum ) = jtile
435     ilocmem ( ichknum ) = itilepos
436     jlocmem ( ichknum ) = jtilepos
437     klocmem ( ichknum ) = layer
438 heimbach 1.8 iobcsmem ( ichknum ) = obcspos
439 heimbach 1.7 icompmem ( ichknum ) = icomp
440     ichkmem ( ichknum ) = ichknum
441     itestmem ( ichknum ) = itest
442     ierrmem ( ichknum ) = ierr
443 heimbach 1.19 icglomem ( ichknum ) = icglo
444    
445 heimbach 1.23 WRITE(standardmessageunit,'(A)')
446     & 'grad-res -------------------------------'
447     WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
448     & ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
449 heimbach 1.19 & layer,itile,jtile,obcspos,
450 heimbach 1.7 & fcref, fcpertplus, fcpertminus
451     #ifdef ALLOW_TANGENTLINEAR_RUN
452 heimbach 1.23 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
453     & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
454 heimbach 1.19 & icompmem(ichknum),itestmem(ichknum),
455     & bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
456 heimbach 1.7 & ftlxxmemo, gfd, ratio_ftl
457 heimbach 1.12 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
458     & 'precision_grdchk_result TLM ', fcref, ftlxxmemo
459     CALL PRINT_MESSAGE
460     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
461 heimbach 1.24 c
462     WRITE(msgBuf,'(A34,1PE24.14)')
463 jmc 1.25 & ' TLM precision_derivative_cost =', fcref
464 heimbach 1.24 CALL PRINT_MESSAGE
465     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
466     WRITE(msgBuf,'(A34,1PE24.14)')
467 jmc 1.25 & ' TLM precision_derivative_grad =', ftlxxmemo
468 heimbach 1.24 CALL PRINT_MESSAGE
469     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
470 heimbach 1.7 #else
471 heimbach 1.23 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
472     & ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
473 heimbach 1.19 & icompmem(ichknum),itestmem(ichknum),
474     & bimem(ichknum),bjmem(ichknum),obcspos,
475 heimbach 1.7 & adxxmemo, gfd, ratio_ad
476 jmc 1.26 c WRITE(msgBuf,'(A34,2(1PE24.14,X))')
477     c & 'precision_grdchk_result ADM ', fcref, adxxmemo
478     c CALL PRINT_MESSAGE
479     c & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
480 heimbach 1.24 WRITE(msgBuf,'(A34,1PE24.14)')
481 jmc 1.25 & ' ADM precision_derivative_cost =', fcref
482 heimbach 1.24 CALL PRINT_MESSAGE
483     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
484     WRITE(msgBuf,'(A34,1PE24.14)')
485 jmc 1.25 & ' ADM precision_derivative_grad =', adxxmemo
486 heimbach 1.24 CALL PRINT_MESSAGE
487     & (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
488 heimbach 1.7 #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 heimbach 1.2 endif
504    
505 heimbach 1.7 c-- end of do icomp = ...
506     enddo
507 heimbach 1.2
508 heimbach 1.7 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 heimbach 1.2
517 heimbach 1.7 c-- Everyone has to wait for the component to be reset.
518     _BARRIER
519 heimbach 1.2
520     c-- Print the results of the gradient check.
521     call grdchk_print( ichknum, ierr_grdchk, mythid )
522    
523 heimbach 1.11 #endif /* ALLOW_GRDCHK */
524 heimbach 1.2
525     end

  ViewVC Help
Powered by ViewVC 1.1.22