/[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.3.6.4 by heimbach, Fri Mar 7 04:01:59 2003 UTC revision 1.33 by jmc, Mon Sep 26 12:56:52 2011 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "GRDCHK_OPTIONS.h"
5    
6  CBOI  CBOI
7  C  C
# Line 30  c         |-- grdchk_loc      - determin Line 31  c         |-- grdchk_loc      - determin
31  c         |  c         |
32  c         |-- grdchk_getxx    - get control vector component from file  c         |-- grdchk_getxx    - get control vector component from file
33  c         |                     perturb it and write back to file  c         |                     perturb it and write back to file
34  c         |-- grdchk_getadxx  - get gradient component calculated  c         |-- grdchk_getadxx  - get gradient component calculated
35  c         |                     via adjoint  c         |                     via adjoint
36  c         |-- the_main_loop   - forward run and cost evaluation  c         |-- the_main_loop   - forward run and cost evaluation
37  c         |                     with perturbed control vector element  c         |                     with perturbed control vector element
# Line 63  c              heimbach@mit.edu 24-Feb-2 Line 64  c              heimbach@mit.edu 24-Feb-2
64  c              - added tangent linear gradient checks  c              - added tangent linear gradient checks
65  c              - fixes for multiproc. gradient checks  c              - fixes for multiproc. gradient checks
66  c              - added more control variables  c              - added more control variables
67  c      c
68  c     ==================================================================  c     ==================================================================
69  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
70  c     ==================================================================  c     ==================================================================
# Line 78  c     == global variables == Line 79  c     == global variables ==
79  #include "PARAMS.h"  #include "PARAMS.h"
80  #include "grdchk.h"  #include "grdchk.h"
81  #include "cost.h"  #include "cost.h"
82    #include "ctrl.h"
83    #ifdef ALLOW_TANGENTLINEAR_RUN
84    #include "g_cost.h"
85    #endif
86    
87  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
88  c     == routine arguments ==  c     == routine arguments ==
89        integer mythid        integer mythid
90    
91  #ifdef ALLOW_GRADIENT_CHECK  #ifdef ALLOW_GRDCHK
92  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
93  c     == local variables ==  c     == local variables ==
94        integer myiter        integer myiter
95        _RL     mytime        _RL     mytime
96        integer bi, itlo, ithi        integer bi, itlo, ithi
97        integer bj, jtlo, jthi        integer bj, jtlo, jthi
98        integer i,  imin, imax        integer i,  imin, imax
# Line 100  c     == local variables == Line 105  c     == local variables ==
105        integer jtile        integer jtile
106        integer itile        integer itile
107        integer layer        integer layer
108          integer obcspos
109        integer itilepos        integer itilepos
110        integer jtilepos        integer jtilepos
111          integer icglo
112        integer itest        integer itest
113        integer ierr        integer ierr
114        integer ierr_grdchk        integer ierr_grdchk
# Line 121  c     == local variables == Line 128  c     == local variables ==
128        _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)
129        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130    
131  #ifdef ALLOW_TANGENTLINEAR_RUN        CHARACTER*(MAX_LEN_MBUF) msgBuf
       _RL     g_fc  
       common /g_cost_r/ g_fc  
 #endif  
132    
133  c     == end of interface ==  c     == end of interface ==
134  CEOP  CEOP
# Line 139  c--   Set the loop ranges. Line 143  c--   Set the loop ranges.
143        imin = 1        imin = 1
144        imax = snx        imax = snx
145    
146          print *, 'ph-check entering grdchk_main '
147    
148  c--   initialise variables  c--   initialise variables
149        call grdchk_init( mythid )        call grdchk_init( mythid )
150    
151  c--   Compute the adjoint models' gradients.  c--   Compute the adjoint model gradients.
152  c--   Compute the unperturbed cost function.  c--   Compute the unperturbed cost function.
153  cph   Gradient via adjoint has already been computed,  cph   Gradient via adjoint has already been computed,
154  cph   and so has unperturbed cost function,  cph   and so has unperturbed cost function,
155  cph   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
156    
157        fcref = fc        ierr      = 0
158        ierr_grdchk = 0        ierr_grdchk = 0
159          adxxmemo  = 0.
160          ftlxxmemo = 0.
161    #ifdef ALLOW_ADMTLM
162          fcref = objf_state_final(idep,jdep,1,1,1)
163    #else
164          fcref = fc
165    #endif
166    
167          print *, 'ph-check fcref = ', fcref
168    
169        do bj = jtlo, jthi        do bj = jtlo, jthi
170           do bi = itlo, ithi           do bi = itlo, ithi
# Line 171  cph   assuming all xx_ fields are initia Line 186  cph   assuming all xx_ fields are initia
186           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
187        end if        end if
188    
189        print *, 'ph-grd 3 -------------------------------'        WRITE(standardmessageunit,'(A)')
190        print ('(2a)'),       &    'grad-res -------------------------------'
191       &     ' ph-grd 3  proc    #    i    j    k            fc ref',        WRITE(standardmessageunit,'(2a)')
192       &     '        fc + eps        fc - eps'       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
193         &    '               fc ref           fc + eps           fc - eps'
194  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
195        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
196       &     ' ph-grd 3  proc    #    i    j    k          tlm grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
197       &     '         fd grad      1 - fd/tlm'       &    '             tlm grad            fd grad         1 - fd/tlm'
198  #else  #else
199        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
200       &     ' ph-grd 3  proc    #    i    j    k          adj grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
201       &     '         fd grad      1 - fd/adj'       &    '             adj grad            fd grad         1 - fd/adj'
202  #endif  #endif
203    
204  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
205  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
206  c--   gradient checks.  c--   gradient checks.
207    
208          if ( nbeg .EQ. 0 )
209         &     call grdchk_get_position( mythid )
210    
211        do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
212    
213           ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
214             adxxmemo  = 0.
215    
216    cph(
217    cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',
218    cph-print     &        nbeg, icomp, ichknum
219    cph)
220           if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
221    
222  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
223              if ( myProcId .EQ. grdchkwhichproc ) then              if ( myProcId .EQ. grdchkwhichproc ) then
224                 call grdchk_loc( icomp, ichknum,                 call grdchk_loc( icomp, ichknum,
225       &              icvrec, itile, jtile, layer,       &              icvrec, itile, jtile, layer, obcspos,
226       &              itilepos, jtilepos, itest, ierr,       &              itilepos, jtilepos, icglo, itest, ierr,
227       &              mythid )       &              mythid )
228    cph(
229    cph-print               print *, 'ph-grd ----- back from loc -----',
230    cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
231    cph)
232              endif              endif
233              _BARRIER              _BARRIER
234                  
235  c******************************************************  c******************************************************
236  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
237  c******************************************************  c******************************************************
# Line 216  c--   get gradient component calculated Line 244  c--   get gradient component calculated
244       &              itilepos, jtilepos,       &              itilepos, jtilepos,
245       &              adxxmemo, mythid )       &              adxxmemo, mythid )
246              endif              endif
247              _BARRIER  C--   Add a global-sum call so that all proc will get the adjoint gradient
248                _GLOBAL_SUM_RL( adxxmemo, myThid )
249    c           _BARRIER
250    
251  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
252  c******************************************************  c******************************************************
# Line 240  c-- Line 270  c--
270  c--   2. perform tangent linear run  c--   2. perform tangent linear run
271              mytime = starttime              mytime = starttime
272              myiter = niter0              myiter = niter0
273    #ifdef ALLOW_ADMTLM
274                do k=1,4*Nr+1
275                 do j=1,sny
276                  do i=1,snx
277                   g_objf_state_final(i,j,1,1,k) = 0.
278                  enddo
279                 enddo
280                enddo
281    #else
282              g_fc = 0.              g_fc = 0.
283    #endif
284    
285              call g_the_main_loop( mytime, myiter, mythid )              call g_the_main_loop( mytime, myiter, mythid )
             ftlxxmemo = g_fc  
286              _BARRIER              _BARRIER
287    #ifdef ALLOW_ADMTLM
288                ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
289    #else
290                ftlxxmemo = g_fc
291    #endif
292    
293  c--  c--
294  c--   3. reset control vector  c--   3. reset control vector
295              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 275  c--   positive perturbation Line 321  c--   positive perturbation
321       &              mythid )       &              mythid )
322              endif              endif
323              _BARRIER              _BARRIER
324                        
325  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
326              mytime = starttime              mytime = starttime
327              myiter = niter0              myiter = niter0
328              call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
329    #ifdef ALLOW_ADMTLM
330                fcpertplus = objf_state_final(idep,jdep,1,1,1)
331    #else
332              fcpertplus = fc              fcpertplus = fc
333    #endif
334                print *, 'ph-check fcpertplus = ', fcpertplus
335              _BARRIER              _BARRIER
336                      
337  c--   Reset control vector.  c--   Reset control vector.
338              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
339       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
# Line 294  c--   Reset control vector. Line 345  c--   Reset control vector.
345              _BARRIER              _BARRIER
346    
347              fcpertminus = fcref              fcpertminus = fcref
348                print *, 'ph-check fcpertminus = ', fcpertminus
349    
350              if ( useCentralDiff ) then              if ( useCentralDiff ) then
351    
# Line 310  c--   repeat the proceedure for a negati Line 362  c--   repeat the proceedure for a negati
362       &                 mythid )       &                 mythid )
363                 endif                 endif
364                 _BARRIER                 _BARRIER
365                        
366  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
367                 mytime = starttime                 mytime = starttime
368                 myiter = niter0                 myiter = niter0
369                 call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
370                 _BARRIER                 _BARRIER
371    #ifdef ALLOW_ADMTLM
372                   fcpertminus = objf_state_final(idep,jdep,1,1,1)
373    #else
374                 fcpertminus = fc                 fcpertminus = fc
375                      #endif
376    
377  c--   Reset control vector.  c--   Reset control vector.
378                 if ( myProcId .EQ. grdchkwhichproc .AND.                 if ( myProcId .EQ. grdchkwhichproc .AND.
379       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
# Line 335  c*************************************** Line 391  c***************************************
391  c--   (D): calculate relative differences between gradients  c--   (D): calculate relative differences between gradients
392  c******************************************************  c******************************************************
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              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
414       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
415    
                if ( grdchk_eps .eq. 0. ) then  
                   gfd = (fcpertplus-fcpertminus)  
                else  
                   gfd = (fcpertplus-fcpertminus)  
      &                 /(grdchk_epsfac*grdchk_eps)  
                endif  
                       
                if ( adxxmemo .eq. 0. ) then  
                   ratio_ad = abs( adxxmemo - gfd )  
                else  
                   ratio_ad = 1. - gfd/adxxmemo  
                endif  
   
                if ( ftlxxmemo .eq. 0. ) then  
                   ratio_ftl = abs( ftlxxmemo - gfd )  
                else  
                   ratio_ftl = 1. - gfd/ftlxxmemo  
                endif  
                     
416                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
417       &              = gfd       &              = gfd
418                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
# Line 381  c*************************************** Line 437  c***************************************
437                 ilocmem   ( ichknum ) = itilepos                 ilocmem   ( ichknum ) = itilepos
438                 jlocmem   ( ichknum ) = jtilepos                 jlocmem   ( ichknum ) = jtilepos
439                 klocmem   ( ichknum ) = layer                 klocmem   ( ichknum ) = layer
440                   iobcsmem  ( ichknum ) = obcspos
441                 icompmem  ( ichknum ) = icomp                 icompmem  ( ichknum ) = icomp
442                 ichkmem   ( ichknum ) = ichknum                 ichkmem   ( ichknum ) = ichknum
443                 itestmem  ( ichknum ) = itest                 itestmem  ( ichknum ) = itest
444                 ierrmem   ( ichknum ) = ierr                 ierrmem   ( ichknum ) = ierr
445                                     icglomem  ( ichknum ) = icglo
446                 print *, 'ph-grd 3 -------------------------------'  
447                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(A)')
448       &              myprocid,ichknum,itilepos,jtilepos,layer,       &             'grad-res -------------------------------'
449                   WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
450         &              ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
451         &              layer,itile,jtile,obcspos,
452       &              fcref, fcpertplus, fcpertminus       &              fcref, fcpertplus, fcpertminus
453  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
454                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
455       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
456       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
457         &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
458       &              ftlxxmemo, gfd, ratio_ftl       &              ftlxxmemo, gfd, ratio_ftl
459  #else  #else
460                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
461       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
462       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
463         &              bimem(ichknum),bjmem(ichknum),obcspos,
464       &              adxxmemo, gfd, ratio_ad       &              adxxmemo, gfd, ratio_ad
465  #endif  #endif
   
466              endif              endif
467    #ifdef ALLOW_TANGENTLINEAR_RUN
468                WRITE(msgBuf,'(A34,1PE24.14)')
469         &              ' TLM  precision_derivative_cost =', fcref
470                CALL PRINT_MESSAGE
471         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
472                WRITE(msgBuf,'(A34,1PE24.14)')
473         &              ' TLM  precision_derivative_grad =', ftlxxmemo
474                CALL PRINT_MESSAGE
475         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
476    #else
477                WRITE(msgBuf,'(A30,1PE22.14)')
478         &              ' ADM  ref_cost_function      =', fcref
479                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
480         &                          SQUEEZE_RIGHT, myThid )
481                WRITE(msgBuf,'(A30,1PE22.14)')
482         &              ' ADM  adjoint_gradient       =', adxxmemo
483                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
484         &                          SQUEEZE_RIGHT, myThid )
485                WRITE(msgBuf,'(A30,1PE22.14)')
486         &              ' ADM  finite-diff_grad       =', gfd
487                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
488         &                          SQUEEZE_RIGHT, myThid )
489    #endif
490    
491              print *, 'ph-grd  ierr ---------------------------'              print *, 'ph-grd  ierr ---------------------------'
492              print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,              print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
# Line 414  c-- else of if ( ichknum ... Line 498  c-- else of if ( ichknum ...
498           else           else
499              ierr_grdchk = -1              ierr_grdchk = -1
500    
501  c-- end of if ( ichknum ...  c-- end of if ( ichknum ...
502           endif           endif
503    
504  c-- end of do icomp = ...  c-- end of do icomp = ...
505        enddo        enddo
506    
507        if ( myProcId .EQ. grdchkwhichproc ) then        if ( myProcId .EQ. grdchkwhichproc ) then
508           CALL WRITE_REC_XYZ_RL(           CALL WRITE_REC_XYZ_RL(
509       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
510           CALL WRITE_REC_XYZ_RL(           CALL WRITE_REC_XYZ_RL(
511       &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)       &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
512           CALL WRITE_REC_XYZ_RL(           CALL WRITE_REC_XYZ_RL(
513       &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)       &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
514        endif        endif
515    
# Line 435  c--   Everyone has to wait for the compo Line 519  c--   Everyone has to wait for the compo
519  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
520        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
521    
522  #endif /* ALLOW_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
523    
524          return
525        end        end
Legend:
Removed from v.1.3.6.4  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22