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

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

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

revision 1.2 by heimbach, Fri Jul 13 14:50:46 2001 UTC revision 1.9 by heimbach, Thu Oct 2 21:30:22 2003 UTC
# Line 2  C $Header$ Line 2  C $Header$
2    
3  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5    CBOI
6    C
7    C !TITLE: GRADIENT CHECK
8    C !AUTHORS: mitgcm developers ( support@mitgcm.org )
9    C !AFFILIATION: Massachussetts Institute of Technology
10    C !DATE:
11    C !INTRODUCTION: gradient check package
12    c \bv
13    c Compare the gradients calculated by the adjoint model with
14    c finite difference approximations.
15    c
16    C     !CALLING SEQUENCE:
17    c
18    c the_model_main
19    c |
20    c |-- ctrl_unpack
21    c |-- adthe_main_loop          - unperturbed cost function and
22    c |-- ctrl_pack                  adjoint gradient are computed here
23    c |
24    c |-- grdchk_main
25    c     |
26    c     |-- grdchk_init
27    c     |-- do icomp=...        - loop over control vector elements
28    c         |
29    c         |-- grdchk_loc      - determine location of icomp on grid
30    c         |
31    c         |-- grdchk_getxx    - get control vector component from file
32    c         |                     perturb it and write back to file
33    c         |-- grdchk_getadxx  - get gradient component calculated
34    c         |                     via adjoint
35    c         |-- the_main_loop   - forward run and cost evaluation
36    c         |                     with perturbed control vector element
37    c         |-- calculate ratio of adj. vs. finite difference gradient
38    c         |
39    c         |-- grdchk_setxx    - Reset control vector element
40    c         |
41    c         |-- grdchk_print    - print results
42    c \ev
43    CEOI
44    
45    CBOP
46    C     !ROUTINE: grdchk_main
47    C     !INTERFACE:
48          subroutine grdchk_main( mythid )
49    
50        subroutine grdchk_main(  C     !DESCRIPTION: \bv
      I                        mythid  
      &                      )  
   
51  c     ==================================================================  c     ==================================================================
52  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
53  c     ==================================================================  c     ==================================================================
 c  
54  c     o Compare the gradients calculated by the adjoint model with  c     o Compare the gradients calculated by the adjoint model with
55  c       finite difference approximations.  c       finite difference approximations.
 c  
56  c     started: Christian Eckert eckert@mit.edu 24-Feb-2000  c     started: Christian Eckert eckert@mit.edu 24-Feb-2000
57  c     continued: heimbach@mit.edu: 13-Jun-2001  c     continued&finished: heimbach@mit.edu: 13-Jun-2001
58  c  c     changed: mlosch@ocean.mit.edu: 09-May-2002
59    c              - added centered difference vs. 1-sided difference option
60    c              - improved output format for readability
61    c              - added control variable hFacC
62    c              heimbach@mit.edu 24-Feb-2003
63    c              - added tangent linear gradient checks
64    c              - fixes for multiproc. gradient checks
65    c              - added more control variables
66    c    
67  c     ==================================================================  c     ==================================================================
68  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
69  c     ==================================================================  c     ==================================================================
70    C     \ev
71    
72    C     !USES:
73        implicit none        implicit none
74    
75  c     == global variables ==  c     == global variables ==
   
76  #include "SIZE.h"  #include "SIZE.h"
77  #include "EEPARAMS.h"  #include "EEPARAMS.h"
78  #include "PARAMS.h"  #include "PARAMS.h"
   
79  #include "grdchk.h"  #include "grdchk.h"
80  #include "cost.h"  #include "cost.h"
81    #ifdef ALLOW_TANGENTLINEAR_RUN
82    #include "g_cost.h"
83    #endif
84    
85    C     !INPUT/OUTPUT PARAMETERS:
86  c     == routine arguments ==  c     == routine arguments ==
   
87        integer mythid        integer mythid
88    
89  #ifdef ALLOW_GRADIENT_CHECK  #ifdef ALLOW_GRADIENT_CHECK
90    C     !LOCAL VARIABLES:
91  c     == local variables ==  c     == local variables ==
92        integer myiter        integer myiter
93        _RL     mytime        _RL     mytime
# Line 46  c     == local variables == Line 97  c     == local variables ==
97        integer j,  jmin, jmax        integer j,  jmin, jmax
98        integer k        integer k
99    
       integer jprocs  
       integer proc_grdchk  
100        integer icomp        integer icomp
101        integer ichknum        integer ichknum
102        integer icvrec        integer icvrec
103        integer jtile        integer jtile
104        integer itile        integer itile
105        integer layer        integer layer
106          integer obcspos
107        integer itilepos        integer itilepos
108        integer jtilepos        integer jtilepos
109        integer itest        integer itest
# Line 61  c     == local variables == Line 111  c     == local variables ==
111        integer ierr_grdchk        integer ierr_grdchk
112        _RL     gfd        _RL     gfd
113        _RL     fcref        _RL     fcref
114        _RL     fcpert        _RL     fcpertplus, fcpertminus
115        _RL     ratio        _RL     ratio_ad
116          _RL     ratio_ftl
117        _RL     xxmemo_ref        _RL     xxmemo_ref
118        _RL     xxmemo_pert        _RL     xxmemo_pert
119        _RL     adxxmemo        _RL     adxxmemo
120          _RL     ftlxxmemo
121          _RL     localEps
122          _RL     grdchk_epsfac
123    
 cph(  
124        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
125        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
126  cph)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
127    
128  c     == end of interface ==  c     == end of interface ==
129    CEOP
130    
131  c--   Set the loop ranges.  c--   Set the loop ranges.
132        jtlo = mybylo(mythid)        jtlo = mybylo(mythid)
# Line 93  cph   Gradient via adjoint has already b Line 147  cph   Gradient via adjoint has already b
147  cph   and so has unperturbed cost function,  cph   and so has unperturbed cost function,
148  cph   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
149    
       fcref = fc  
150        ierr_grdchk = 0        ierr_grdchk = 0
151    cphadmtlm(
152          fcref = fc
153    cphadmtlm      fcref = objf_state_final(45,4,1,1)
154    cphadmtlm)
155    
156          print *, 'ph-check fcref = ', fcref
157    
 cph(  
158        do bj = jtlo, jthi        do bj = jtlo, jthi
159           do bi = itlo, ithi           do bi = itlo, ithi
160              do k = 1, nr              do k = 1, nr
# Line 104  cph( Line 162  cph(
162                    do i = imin, imax                    do i = imin, imax
163                       tmpplot1(i,j,k,bi,bj) = 0. _d 0                       tmpplot1(i,j,k,bi,bj) = 0. _d 0
164                       tmpplot2(i,j,k,bi,bj) = 0. _d 0                       tmpplot2(i,j,k,bi,bj) = 0. _d 0
165                         tmpplot3(i,j,k,bi,bj) = 0. _d 0
166                    end do                    end do
167                 end do                 end do
168              end do              end do
169           end do           end do
170        end do        end do
171  cph)  
172          if ( useCentralDiff ) then
173             grdchk_epsfac = 2. _d 0
174          else
175             grdchk_epsfac = 1. _d 0
176          end if
177    
178          print *, 'ph-grd 3 -------------------------------'
179          print ('(2a)'),
180         &     ' ph-grd 3  proc    #    i    j    k            fc ref',
181         &     '        fc + eps        fc - eps'
182    #ifdef ALLOW_TANGENTLINEAR_RUN
183          print ('(2a)'),
184         &     ' ph-grd 3  proc    #    i    j    k          tlm grad',
185         &     '         fd grad      1 - fd/tlm'
186    #else
187          print ('(2a)'),
188         &     ' ph-grd 3  proc    #    i    j    k          adj grad',
189         &     '         fd grad      1 - fd/adj'
190    #endif
191    
192  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
193  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
194  c--   gradient checks.  c--   gradient checks.
       do jprocs = 1,numberOfProcs  
          proc_grdchk = jprocs - 1  
   
          if ( myProcId .eq. proc_grdchk ) then  
195    
196              do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
197    
198                 ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
199    
200                 if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
201    
202  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
203                    call grdchk_loc( icomp, ichknum,              if ( myProcId .EQ. grdchkwhichproc ) then
204       &                 icvrec, itile, jtile, layer,                 call grdchk_loc( icomp, ichknum,
205       &                 itilepos, jtilepos, itest, ierr,       &              icvrec, itile, jtile, layer, obcspos,
206       &                 mythid )       &              itilepos, jtilepos, itest, ierr,
207         &              mythid )
208                endif
209                _BARRIER
210                                
211                    if ( ierr .eq. 0 ) then  c******************************************************
212    c--   (A): get gradient component calculated via adjoint
213    c******************************************************
214    
215    c--   get gradient component calculated via adjoint
216                if ( myProcId .EQ. grdchkwhichproc .AND.
217         &           ierr .EQ. 0 ) then
218                   call grdchk_getadxx( icvrec,
219         &              itile, jtile, layer,
220         &              itilepos, jtilepos,
221         &              adxxmemo, mythid )
222                endif
223                _BARRIER
224    
225    #ifdef ALLOW_TANGENTLINEAR_RUN
226    c******************************************************
227    c--   (B): Get gradient component g_fc from tangent linear run:
228    c******************************************************
229    c--
230    c--   1. perturb control vector component: xx(i)=1.
231    
232                if ( myProcId .EQ. grdchkwhichproc .AND.
233         &           ierr .EQ. 0 ) then
234                   localEps = 1. _d 0
235                   call grdchk_getxx( icvrec, TANGENT_SIMULATION,
236         &              itile, jtile, layer,
237         &              itilepos, jtilepos,
238         &              xxmemo_ref, xxmemo_pert, localEps,
239         &              mythid )
240                endif
241                _BARRIER
242    
243    c--
244    c--   2. perform tangent linear run
245                mytime = starttime
246                myiter = niter0
247    cphadmtlm(
248                g_fc = 0.
249    cphadmtlm            do j=1,sny
250    cphadmtlm               do i=1,snx
251    cphadmtlm                  g_objf_state_final(i,j,1,1) = 0.
252    cphadmtlm               enddo
253    cphadmtlm            enddo
254    cphadmtlm)
255                call g_the_main_loop( mytime, myiter, mythid )
256    cphadmtlm(
257                ftlxxmemo = g_fc
258    cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1)
259    cphadmtlm)
260                _BARRIER
261    c--
262    c--   3. reset control vector
263                if ( myProcId .EQ. grdchkwhichproc .AND.
264         &           ierr .EQ. 0 ) then
265                   call grdchk_setxx( icvrec, TANGENT_SIMULATION,
266         &              itile, jtile, layer,
267         &              itilepos, jtilepos,
268         &              xxmemo_ref, mythid )
269                endif
270                _BARRIER
271    
272    #endif /* ALLOW_TANGENTLINEAR_RUN */
273    
274    
275    c******************************************************
276    c--   (C): Get gradient via finite difference perturbation
277    c******************************************************
278    
279  c--   get control vector component from file  c--   get control vector component from file
280  c--   perturb it and write back to file  c--   perturb it and write back to file
281                       call grdchk_getxx( icvrec,  c--   positive perturbation
282       &                    itile, jtile, layer,              localEps = abs(grdchk_eps)
283       &                    itilepos, jtilepos,              if ( myProcId .EQ. grdchkwhichproc .AND.
284       &                    xxmemo_ref, xxmemo_pert, mythid )       &           ierr .EQ. 0 ) then
285                       _BARRIER                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
286         &              itile, jtile, layer,
287         &              itilepos, jtilepos,
288         &              xxmemo_ref, xxmemo_pert, localEps,
289         &              mythid )
290                endif
291                _BARRIER
292                                            
 c--   get gradient component calculated via adjoint  
                      call grdchk_getadxx( icvrec,  
      &                    itile, jtile, layer,  
      &                    itilepos, jtilepos,  
      &                    adxxmemo, mythid )  
                      _BARRIER  
   
293  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
294                       mytime = starttime              mytime = starttime
295                       myiter = niter0              myiter = niter0
296                       call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
297                       fcpert = fc  cphadmtlm(
298                fcpertplus = fc
299    cphadmtlm            fcpertplus = objf_state_final(45,4,1,1)
300    cphadmtlm)
301                _BARRIER
302                                        
303                       if ( grdchk_eps .eq. 0. ) then  c--   Reset control vector.
304                          gfd = (fcpert-fcref)              if ( myProcId .EQ. grdchkwhichproc .AND.
305                       else       &           ierr .EQ. 0 ) then
306                          gfd = (fcpert-fcref)/grdchk_eps                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
307                       endif       &              itile, jtile, layer,
308         &              itilepos, jtilepos,
309         &              xxmemo_ref, mythid )
310                endif
311                _BARRIER
312    
313                fcpertminus = fcref
314    
315                if ( useCentralDiff ) then
316    
317    c--   get control vector component from file
318    c--   perturb it and write back to file
319    c--   repeat the proceedure for a negative perturbation
320                   if ( myProcId .EQ. grdchkwhichproc .AND.
321         &           ierr .EQ. 0 ) then
322                      localEps = - abs(grdchk_eps)
323                      call grdchk_getxx( icvrec, FORWARD_SIMULATION,
324         &                 itile, jtile, layer,
325         &                 itilepos, jtilepos,
326         &                 xxmemo_ref, xxmemo_pert, localEps,
327         &                 mythid )
328                   endif
329                   _BARRIER
330                                            
331                       if ( gfd .eq. 0. ) then  c--   forward run with perturbed control vector
332                          ratio = abs( adxxmemo - gfd )                 mytime = starttime
333                       else                 myiter = niter0
334                          ratio = 1. - adxxmemo/gfd                 call the_main_loop( mytime, myiter, mythid )
335                       endif                 _BARRIER
336                   fcpertminus = fc
337                                        
338  c--   Reset control vector.  c--   Reset control vector.
339                       call grdchk_setxx( icvrec,                 if ( myProcId .EQ. grdchkwhichproc .AND.
340       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
341       &                    itilepos, jtilepos,                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,
342       &                    xxmemo_ref, mythid )       &                 itile, jtile, layer,
343                       _BARRIER       &                 itilepos, jtilepos,
344         &                 xxmemo_ref, mythid )
345  cph(                 endif
346                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =                 _BARRIER
      &                    gfd  
                      tmpplot2(itilepos,jtilepos,layer,itile,jtile) =  
      &                    ratio  
 cph)  
347    
348    c-- end of if useCentralDiff ...
349                end if
350    
351    c******************************************************
352    c--   (D): calculate relative differences between gradients
353    c******************************************************
354    
355                if ( myProcId .EQ. grdchkwhichproc .AND.
356         &           ierr .EQ. 0 ) then
357    
358                   if ( grdchk_eps .eq. 0. ) then
359                      gfd = (fcpertplus-fcpertminus)
360                   else
361                      gfd = (fcpertplus-fcpertminus)
362         &                 /(grdchk_epsfac*grdchk_eps)
363                   endif
364                                            
365                       fcpmem    ( ichknum ) = fcpert                 if ( adxxmemo .eq. 0. ) then
366                       xxmemref  ( ichknum ) = xxmemo_ref                    ratio_ad = abs( adxxmemo - gfd )
                      xxmempert ( ichknum ) = xxmemo_pert  
                      gfdmem    ( ichknum ) = gfd  
                      adxxmem   ( ichknum ) = adxxmemo  
                      ratiomem  ( ichknum ) = ratio  
   
                      irecmem   ( ichknum ) = icvrec  
                      bimem     ( ichknum ) = itile  
                      bjmem     ( ichknum ) = jtile  
                      ilocmem   ( ichknum ) = itilepos  
                      jlocmem   ( ichknum ) = jtilepos  
                      klocmem   ( ichknum ) = layer  
                      icompmem  ( ichknum ) = icomp  
                      ichkmem   ( ichknum ) = ichknum  
                      itestmem  ( ichknum ) = itest  
                      ierrmem   ( ichknum ) = ierr  
                     
 cph(  
                      print *, 'ph-grd 3 -------------------------------'  
                      print *, 'ph-grd 3 ',  
      &                    ichknum,itilepos,jtilepos,layer,  
      &                    fcref, fcpert  
                      print *, 'ph-grd 3 ',  
      &                    ichknum,ichkmem(ichknum),  
      &                    icompmem(ichknum),itestmem(ichknum),  
      &                    adxxmemo, gfd, ratio  
 cph)  
                   else  
 c  
                   endif  
367                 else                 else
368                    ierr_grdchk = -1                    ratio_ad = 1. - gfd/adxxmemo
369                 endif                 endif
               
             enddo  
          endif  
370    
371  c--   Everyone has to wait for the component to be reset.                 if ( ftlxxmemo .eq. 0. ) then
372           _BARRIER                    ratio_ftl = abs( ftlxxmemo - gfd )
373                   else
374                      ratio_ftl = 1. - gfd/ftlxxmemo
375                   endif
376                      
377                   tmpplot1(itilepos,jtilepos,layer,itile,jtile)
378         &              = gfd
379                   tmpplot2(itilepos,jtilepos,layer,itile,jtile)
380         &              = ratio_ad
381                   tmpplot3(itilepos,jtilepos,layer,itile,jtile)
382         &              = ratio_ftl
383    
384                   fcrmem      ( ichknum ) = fcref
385                   fcppmem     ( ichknum ) = fcpertplus
386                   fcpmmem     ( ichknum ) = fcpertminus
387                   xxmemref    ( ichknum ) = xxmemo_ref
388                   xxmempert   ( ichknum ) = xxmemo_pert
389                   gfdmem      ( ichknum ) = gfd
390                   adxxmem     ( ichknum ) = adxxmemo
391                   ftlxxmem    ( ichknum ) = ftlxxmemo
392                   ratioadmem  ( ichknum ) = ratio_ad
393                   ratioftlmem ( ichknum ) = ratio_ftl
394    
395                   irecmem   ( ichknum ) = icvrec
396                   bimem     ( ichknum ) = itile
397                   bjmem     ( ichknum ) = jtile
398                   ilocmem   ( ichknum ) = itilepos
399                   jlocmem   ( ichknum ) = jtilepos
400                   klocmem   ( ichknum ) = layer
401                   iobcsmem  ( ichknum ) = obcspos
402                   icompmem  ( ichknum ) = icomp
403                   ichkmem   ( ichknum ) = ichknum
404                   itestmem  ( ichknum ) = itest
405                   ierrmem   ( ichknum ) = ierr
406                      
407                   print *, 'ph-grd 3 -------------------------------'
408                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
409         &              myprocid,ichknum,itilepos,jtilepos,layer,
410         &              fcref, fcpertplus, fcpertminus
411    #ifdef ALLOW_TANGENTLINEAR_RUN
412                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
413         &              myprocid,ichknum,ichkmem(ichknum),
414         &              icompmem(ichknum),itestmem(ichknum),
415         &              ftlxxmemo, gfd, ratio_ftl
416    #else
417                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
418         &              myprocid,ichknum,ichkmem(ichknum),
419         &              icompmem(ichknum),itestmem(ichknum),
420         &              adxxmemo, gfd, ratio_ad
421    #endif
422    
423                endif
424    
425                print *, 'ph-grd  ierr ---------------------------'
426                print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
427         &           ', ichknum = ', ichknum
428    
429                _BARRIER
430    
431    c-- else of if ( ichknum ...
432             else
433                ierr_grdchk = -1
434    
435    c-- end of if ( ichknum ...
436             endif
437    
438    c-- end of do icomp = ...
439        enddo        enddo
440    
441  cph(        if ( myProcId .EQ. grdchkwhichproc ) then
442        CALL WRITE_REC_XYZ_RL( 'grd_findiff' , tmpplot1, 1, 0, myThid)           CALL WRITE_REC_XYZ_RL(
443        CALL WRITE_REC_XYZ_RL( 'grd_ratio' , tmpplot2, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
444  cph)           CALL WRITE_REC_XYZ_RL(
445         &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
446             CALL WRITE_REC_XYZ_RL(
447         &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
448          endif
449    
450    c--   Everyone has to wait for the component to be reset.
451          _BARRIER
452    
453  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
454        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22