/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_common/grdchk_main.F
ViewVC logotype

Annotation of /MITgcm_contrib/heimbach/OpenAD/code_common/grdchk_main.F

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


Revision 1.3 - (hide annotations) (download)
Wed Aug 22 22:35:45 2012 UTC (12 years, 11 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
FILE REMOVED
changes roled in main version

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

  ViewVC Help
Powered by ViewVC 1.1.22