/[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.20 by heimbach, Sat Mar 24 02:25:29 2007 UTC revision 1.31 by jmc, Tue May 24 22:40:30 2011 UTC
# Line 1  Line 1 
 C  
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "AD_CONFIG.h"  #include "GRDCHK_OPTIONS.h"
 #include "CPP_OPTIONS.h"  
5    
6  CBOI  CBOI
7  C  C
# Line 33  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 66  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 81  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  #ifdef ALLOW_TANGENTLINEAR_RUN
84  #include "g_cost.h"  #include "g_cost.h"
85  #endif  #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_GRDCHK  #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 149  c--   Set the loop ranges. Line 148  c--   Set the loop ranges.
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          ierr      = 0
158        ierr_grdchk = 0        ierr_grdchk = 0
159  cphadmtlm(        adxxmemo  = 0.
160          ftlxxmemo = 0.
161    #ifdef ALLOW_ADMTLM
162          fcref = objf_state_final(idep,jdep,1,1,1)
163    #else
164        fcref = fc        fcref = fc
165  cphadmtlm      fcref = objf_state_final(45,4,1,1,1)  #endif
 cphadmtlm)  
166    
167        print *, 'ph-check fcref = ', fcref        print *, 'ph-check fcref = ', fcref
168    
# Line 183  cphadmtlm) Line 186  cphadmtlm)
186           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
187        end if        end if
188    
189        print *, 'grad-res -------------------------------'        WRITE(standardmessageunit,'(A)')
190        print '(2a)',       &    'grad-res -------------------------------'
191          WRITE(standardmessageunit,'(2a)')
192       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
193       &    '               fc ref           fc + eps           fc - eps'       &    '               fc ref           fc + eps           fc - eps'
194  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
195        print '(2a)',        WRITE(standardmessageunit,'(2a)')
196       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
197       &    '             tlm grad            fd grad         1 - fd/tlm'       &    '             tlm grad            fd grad         1 - fd/tlm'
198  #else  #else
199        print '(2a)',        WRITE(standardmessageunit,'(2a)')
200       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
201       &    '             adj grad            fd grad         1 - fd/adj'       &    '             adj grad            fd grad         1 - fd/adj'
202  #endif  #endif
# Line 201  c--   Compute the finite difference appr Line 205  c--   Compute the finite difference appr
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 )        if ( nbeg .EQ. 0 )
209       &     call grdchk_get_position( mythid )       &     call grdchk_get_position( mythid )
210    
211        do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
# Line 209  c--   gradient checks. Line 213  c--   gradient checks.
213           ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
214    
215  cph(  cph(
216  cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',  cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',
217  cph-print     &        nbeg, icomp, ichknum  cph-print     &        nbeg, icomp, ichknum
218  cph)  cph)
219           if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
220    
221  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
222              if ( myProcId .EQ. grdchkwhichproc ) then              if ( myProcId .EQ. grdchkwhichproc ) then
223                 call grdchk_loc( icomp, ichknum,                 call grdchk_loc( icomp, ichknum,
224       &              icvrec, itile, jtile, layer, obcspos,       &              icvrec, itile, jtile, layer, obcspos,
225       &              itilepos, jtilepos, icglo, itest, ierr,       &              itilepos, jtilepos, icglo, itest, ierr,
226       &              mythid )       &              mythid )
# Line 226  cph-print     &             icvrec, itil Line 230  cph-print     &             icvrec, itil
230  cph)  cph)
231              endif              endif
232              _BARRIER              _BARRIER
233                  
234  c******************************************************  c******************************************************
235  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
236  c******************************************************  c******************************************************
# Line 263  c-- Line 267  c--
267  c--   2. perform tangent linear run  c--   2. perform tangent linear run
268              mytime = starttime              mytime = starttime
269              myiter = niter0              myiter = niter0
270  cphadmtlm(  #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.              g_fc = 0.
280  cphadmtlm            do j=1,sny  #endif
281  cphadmtlm               do i=1,snx  
 cphadmtlm                  g_objf_state_final(i,j,1,1,1) = 0.  
 cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 0.  
 cphadmtlm               enddo  
 cphadmtlm            enddo  
 cphadmtlm)  
282              call g_the_main_loop( mytime, myiter, mythid )              call g_the_main_loop( mytime, myiter, mythid )
 cphadmtlm(  
             ftlxxmemo = g_fc  
 cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1,1)  
 cphadmtlm)  
283              _BARRIER              _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--  c--
291  c--   3. reset control vector  c--   3. reset control vector
292              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 309  c--   positive perturbation Line 318  c--   positive perturbation
318       &              mythid )       &              mythid )
319              endif              endif
320              _BARRIER              _BARRIER
321                        
322  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
323              mytime = starttime              mytime = starttime
324              myiter = niter0              myiter = niter0
325              call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
326  cphadmtlm(  #ifdef ALLOW_ADMTLM
327                fcpertplus = objf_state_final(idep,jdep,1,1,1)
328    #else
329              fcpertplus = fc              fcpertplus = fc
330  cphadmtlm            fcpertplus = objf_state_final(45,4,1,1,1)  #endif
 cphadmtlm)  
331              print *, 'ph-check fcpertplus = ', fcpertplus              print *, 'ph-check fcpertplus = ', fcpertplus
332              _BARRIER              _BARRIER
333                      
334  c--   Reset control vector.  c--   Reset control vector.
335              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
336       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
# Line 349  c--   repeat the proceedure for a negati Line 359  c--   repeat the proceedure for a negati
359       &                 mythid )       &                 mythid )
360                 endif                 endif
361                 _BARRIER                 _BARRIER
362                        
363  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
364                 mytime = starttime                 mytime = starttime
365                 myiter = niter0                 myiter = niter0
366                 call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
367                 _BARRIER                 _BARRIER
368    #ifdef ALLOW_ADMTLM
369                   fcpertminus = objf_state_final(idep,jdep,1,1,1)
370    #else
371                 fcpertminus = fc                 fcpertminus = fc
372                      #endif
373    
374  c--   Reset control vector.  c--   Reset control vector.
375                 if ( myProcId .EQ. grdchkwhichproc .AND.                 if ( myProcId .EQ. grdchkwhichproc .AND.
376       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
# Line 383  c*************************************** Line 397  c***************************************
397                    gfd = (fcpertplus-fcpertminus)                    gfd = (fcpertplus-fcpertminus)
398       &                 /(grdchk_epsfac*grdchk_eps)       &                 /(grdchk_epsfac*grdchk_eps)
399                 endif                 endif
400                        
401                 if ( adxxmemo .eq. 0. ) then                 if ( adxxmemo .eq. 0. ) then
402                    ratio_ad = abs( adxxmemo - gfd )                    ratio_ad = abs( adxxmemo - gfd )
403                 else                 else
# Line 395  c*************************************** Line 409  c***************************************
409                 else                 else
410                    ratio_ftl = 1. - gfd/ftlxxmemo                    ratio_ftl = 1. - gfd/ftlxxmemo
411                 endif                 endif
412                      
413                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)
414       &              = gfd       &              = gfd
415                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)
# Line 427  c*************************************** Line 441  c***************************************
441                 ierrmem   ( ichknum ) = ierr                 ierrmem   ( ichknum ) = ierr
442                 icglomem  ( ichknum ) = icglo                 icglomem  ( ichknum ) = icglo
443    
444                 print *, 'grad-res -------------------------------'                 WRITE(standardmessageunit,'(A)')
445                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',       &             'grad-res -------------------------------'
446       &              myprocid,ichknum,itilepos,jtilepos,                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
447         &              ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
448       &              layer,itile,jtile,obcspos,       &              layer,itile,jtile,obcspos,
449       &              fcref, fcpertplus, fcpertminus       &              fcref, fcpertplus, fcpertminus
450  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
451                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
452       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
453       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
454       &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),       &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
455       &              ftlxxmemo, gfd, ratio_ftl       &              ftlxxmemo, gfd, ratio_ftl
456                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
457       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
458                 CALL PRINT_MESSAGE                 CALL PRINT_MESSAGE
459       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
460    c
461                   WRITE(msgBuf,'(A34,1PE24.14)')
462         &              ' TLM  precision_derivative_cost =', fcref
463                   CALL PRINT_MESSAGE
464         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
465                   WRITE(msgBuf,'(A34,1PE24.14)')
466         &              ' TLM  precision_derivative_grad =', ftlxxmemo
467                   CALL PRINT_MESSAGE
468         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
469  #else  #else
470                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
471       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
472       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
473       &              bimem(ichknum),bjmem(ichknum),obcspos,       &              bimem(ichknum),bjmem(ichknum),obcspos,
474       &              adxxmemo, gfd, ratio_ad       &              adxxmemo, gfd, ratio_ad
475                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')  c              WRITE(msgBuf,'(A34,2(1PE24.14,X))')
476       &              'precision_grdchk_result ADM ', fcref, adxxmemo  c    &              'precision_grdchk_result ADM ', fcref, adxxmemo
477    c              CALL PRINT_MESSAGE
478    c    &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
479                   WRITE(msgBuf,'(A34,1PE24.14)')
480         &              ' ADM  precision_derivative_cost =', fcref
481                   CALL PRINT_MESSAGE
482         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
483                   WRITE(msgBuf,'(A34,1PE24.14)')
484         &              ' ADM  precision_derivative_grad =', adxxmemo
485                 CALL PRINT_MESSAGE                 CALL PRINT_MESSAGE
486       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
487  #endif  #endif
# Line 466  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 489  c--   Print the results of the gradient Line 521  c--   Print the results of the gradient
521    
522  #endif /* ALLOW_GRDCHK */  #endif /* ALLOW_GRDCHK */
523    
524          return
525        end        end
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22