/[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.6 by heimbach, Mon Sep 16 18:11:58 2002 UTC revision 1.12 by heimbach, Tue Nov 4 20:47:42 2003 UTC
# Line 1  Line 1 
1    C
2  C $Header$  C $Header$
3    C $Name$
4    
5    #include "AD_CONFIG.h"
6  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
7    
8  CBOI  CBOI
# Line 59  c     changed: mlosch@ocean.mit.edu: 09- Line 62  c     changed: mlosch@ocean.mit.edu: 09-
62  c              - added centered difference vs. 1-sided difference option  c              - added centered difference vs. 1-sided difference option
63  c              - improved output format for readability  c              - improved output format for readability
64  c              - added control variable hFacC  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  c      c    
70  c     ==================================================================  c     ==================================================================
71  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
# Line 74  c     == global variables == Line 81  c     == global variables ==
81  #include "PARAMS.h"  #include "PARAMS.h"
82  #include "grdchk.h"  #include "grdchk.h"
83  #include "cost.h"  #include "cost.h"
84    #ifdef ALLOW_TANGENTLINEAR_RUN
85    #include "g_cost.h"
86    #endif
87    
88  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
89  c     == routine arguments ==  c     == routine arguments ==
90        integer mythid        integer mythid
91    
92  #ifdef ALLOW_GRADIENT_CHECK  #ifdef ALLOW_GRDCHK
93  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
94  c     == local variables ==  c     == local variables ==
95        integer myiter        integer myiter
# Line 90  c     == local variables == Line 100  c     == local variables ==
100        integer j,  jmin, jmax        integer j,  jmin, jmax
101        integer k        integer k
102    
       integer jprocs  
       integer proc_grdchk  
103        integer icomp        integer icomp
104        integer ichknum        integer ichknum
105        integer icvrec        integer icvrec
106        integer jtile        integer jtile
107        integer itile        integer itile
108        integer layer        integer layer
109          integer obcspos
110        integer itilepos        integer itilepos
111        integer jtilepos        integer jtilepos
112        integer itest        integer itest
# Line 115  c     == local variables == Line 124  c     == local variables ==
124        _RL     localEps        _RL     localEps
125        _RL     grdchk_epsfac        _RL     grdchk_epsfac
126    
 #ifdef ALLOW_TANGENTLINEAR_RUN  
       _RL     g_fc  
       common /g_cost_r/ g_fc  
 #endif  
   
127        _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)
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          CHARACTER*(MAX_LEN_MBUF) msgBuf
132    
133  c     == end of interface ==  c     == end of interface ==
134  CEOP  CEOP
135    
# Line 142  c--   initialise variables Line 148  c--   initialise variables
148    
149  c--   Compute the adjoint models' gradients.  c--   Compute the adjoint models' gradients.
150  c--   Compute the unperturbed cost function.  c--   Compute the unperturbed cost function.
151  c--   Gradient via adjoint has already been computed,  cph   Gradient via adjoint has already been computed,
152  c--   and so has unperturbed cost function,  cph   and so has unperturbed cost function,
153  c--   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
154    
       fcref = fc  
155        ierr_grdchk = 0        ierr_grdchk = 0
156    cphadmtlm(
157          fcref = fc
158    cphadmtlm      fcref = objf_state_final(45,4,1,1)
159    cphadmtlm)
160    
161          print *, 'ph-check fcref = ', fcref
162    
163        do bj = jtlo, jthi        do bj = jtlo, jthi
164           do bi = itlo, ithi           do bi = itlo, ithi
# Line 169  c--   assuming all xx_ fields are initia Line 180  c--   assuming all xx_ fields are initia
180           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
181        end if        end if
182    
183          print *, 'grad-res -------------------------------'
184          print ('(2a)'),
185         &     ' grad-res  proc    #    i    j    k            fc ref',
186         &     '        fc + eps        fc - eps'
187    #ifdef ALLOW_TANGENTLINEAR_RUN
188          print ('(2a)'),
189         &     ' grad-res  proc    #    i    j    k          tlm grad',
190         &     '         fd grad      1 - fd/tlm'
191    #else
192          print ('(2a)'),
193         &     ' grad-res  proc    #    i    j    k          adj grad',
194         &     '         fd grad      1 - fd/adj'
195    #endif
196    
197  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
198  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
199  c--   gradient checks.  c--   gradient checks.
       do jprocs = 1,numberOfProcs  
          proc_grdchk = jprocs - 1  
   
          if ( myProcId .eq. proc_grdchk ) then  
200    
201              do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
202    
203                 ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
204    
205                 if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
206    
207  c--   Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
208                    call grdchk_loc( icomp, ichknum,              if ( myProcId .EQ. grdchkwhichproc ) then
209       &                 icvrec, itile, jtile, layer,                 call grdchk_loc( icomp, ichknum,
210       &                 itilepos, jtilepos, itest, ierr,       &              icvrec, itile, jtile, layer, obcspos,
211       &                 mythid )       &              itilepos, jtilepos, itest, ierr,
212         &              mythid )
213                endif
214                _BARRIER
215                                
                   if ( ierr .eq. 0 ) then  
   
216  c******************************************************  c******************************************************
217  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
218  c******************************************************  c******************************************************
219                       call grdchk_getadxx( icvrec,  
220       &                    itile, jtile, layer,  c--   get gradient component calculated via adjoint
221       &                    itilepos, jtilepos,              if ( myProcId .EQ. grdchkwhichproc .AND.
222       &                    adxxmemo, mythid )       &           ierr .EQ. 0 ) then
223                       _BARRIER                 call grdchk_getadxx( icvrec,
224         &              itile, jtile, layer,
225         &              itilepos, jtilepos,
226         &              adxxmemo, mythid )
227                endif
228                _BARRIER
229    
230  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
231  c******************************************************  c******************************************************
# Line 207  c*************************************** Line 234  c***************************************
234  c--  c--
235  c--   1. perturb control vector component: xx(i)=1.  c--   1. perturb control vector component: xx(i)=1.
236    
237                       localEps = 1.              if ( myProcId .EQ. grdchkwhichproc .AND.
238                       call grdchk_getxx( icvrec, TANGENT_SIMULATION,       &           ierr .EQ. 0 ) then
239       &                    itile, jtile, layer,                 localEps = 1. _d 0
240       &                    itilepos, jtilepos,                 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
241       &                    xxmemo_ref, xxmemo_pert, localEps,       &              itile, jtile, layer,
242       &                    mythid )       &              itilepos, jtilepos,
243                       _BARRIER       &              xxmemo_ref, xxmemo_pert, localEps,
244         &              mythid )
245                endif
246                _BARRIER
247    
248  c--  c--
249  c--   2. perform tangent linear run  c--   2. perform tangent linear run
250                       mytime = starttime              mytime = starttime
251                       myiter = niter0              myiter = niter0
252                       g_fc = 0.  cphadmtlm(
253                       call g_the_main_loop( mytime, myiter, mythid )              g_fc = 0.
254                       ftlxxmemo = g_fc  cphadmtlm            do j=1,sny
255    cphadmtlm               do i=1,snx
256    cphadmtlm                  g_objf_state_final(i,j,1,1) = 0.
257    cphadmtlm               enddo
258    cphadmtlm            enddo
259    cphadmtlm)
260                call g_the_main_loop( mytime, myiter, mythid )
261    cphadmtlm(
262                ftlxxmemo = g_fc
263    cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1)
264    cphadmtlm)
265                _BARRIER
266  c--  c--
267  c--   3. reset control vector  c--   3. reset control vector
268                       call grdchk_setxx( icvrec, TANGENT_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
269       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
270       &                    itilepos, jtilepos,                 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
271       &                    xxmemo_ref, mythid )       &              itile, jtile, layer,
272                       _BARRIER       &              itilepos, jtilepos,
273  #endif       &              xxmemo_ref, mythid )
274                endif
275                _BARRIER
276    
277    #endif /* ALLOW_TANGENTLINEAR_RUN */
278    
279    
280  c******************************************************  c******************************************************
281  c--   (C): Get gradient via finite difference perturbation  c--   (C): Get gradient via finite difference perturbation
282  c******************************************************  c******************************************************
283    
284  c--   get control vector component from file  c--   get control vector component from file
285  c--   perturb it and write back to file:  c--   perturb it and write back to file
286  c--   positive perturbation  c--   positive perturbation
287                       localEps = abs(grdchk_eps)              localEps = abs(grdchk_eps)
288                       call grdchk_getxx( icvrec, FORWARD_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
289       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
290       &                    itilepos, jtilepos,                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
291       &                    xxmemo_ref, xxmemo_pert, localEps,       &              itile, jtile, layer,
292       &                    mythid )       &              itilepos, jtilepos,
293                       _BARRIER       &              xxmemo_ref, xxmemo_pert, localEps,
294         &              mythid )
295                endif
296                _BARRIER
297                                            
298  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
299                       mytime = starttime              mytime = starttime
300                       myiter = niter0              myiter = niter0
301                       call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
302                       fcpertplus = fc  cphadmtlm(
303                fcpertplus = fc
304    cphadmtlm            fcpertplus = objf_state_final(45,4,1,1)
305    cphadmtlm)
306                _BARRIER
307                                        
308  c--   Reset control vector.  c--   Reset control vector.
309                       call grdchk_setxx( icvrec, FORWARD_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
310       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
311       &                    itilepos, jtilepos,                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
312       &                    xxmemo_ref, mythid )       &              itile, jtile, layer,
313                       _BARRIER       &              itilepos, jtilepos,
314         &              xxmemo_ref, mythid )
315                endif
316                _BARRIER
317    
318                       fcpertminus = fcref              fcpertminus = fcref
319    
320                       if ( useCentralDiff ) then              if ( useCentralDiff ) then
321    
322  c--   get control vector component from file  c--   get control vector component from file
323  c--   perturb it and write back to file:  c--   perturb it and write back to file
324  c--   repeat the proceedure for a negative perturbation  c--   repeat the proceedure for a negative perturbation
325                          localEps = - abs(grdchk_eps)                 if ( myProcId .EQ. grdchkwhichproc .AND.
326                          call grdchk_getxx( icvrec, FORWARD_SIMULATION,       &           ierr .EQ. 0 ) then
327       &                    itile, jtile, layer,                    localEps = - abs(grdchk_eps)
328       &                    itilepos, jtilepos,                    call grdchk_getxx( icvrec, FORWARD_SIMULATION,
329       &                    xxmemo_ref, xxmemo_pert, localEps,       &                 itile, jtile, layer,
330       &                    mythid )       &                 itilepos, jtilepos,
331                          _BARRIER       &                 xxmemo_ref, xxmemo_pert, localEps,
332         &                 mythid )
333                   endif
334                   _BARRIER
335                                            
336  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
337                          mytime = starttime                 mytime = starttime
338                          myiter = niter0                 myiter = niter0
339                          call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
340                          fcpertminus = fc                 _BARRIER
341                   fcpertminus = fc
342                                        
343  c--   Reset control vector.  c--   Reset control vector.
344                          call grdchk_setxx( icvrec, FORWARD_SIMULATION,                 if ( myProcId .EQ. grdchkwhichproc .AND.
345       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
346       &                    itilepos, jtilepos,                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,
347       &                    xxmemo_ref, mythid )       &                 itile, jtile, layer,
348                          _BARRIER       &                 itilepos, jtilepos,
349         &                 xxmemo_ref, mythid )
350                   endif
351                   _BARRIER
352    
353    c-- end of if useCentralDiff ...
354                end if
355    
                      end if  
 c--  
356  c******************************************************  c******************************************************
357  c--   (D): calculate relative differences between gradients  c--   (D): calculate relative differences between gradients
358  c******************************************************  c******************************************************
359    
360                       if ( grdchk_eps .eq. 0. ) then              if ( myProcId .EQ. grdchkwhichproc .AND.
361                          gfd = (fcpertplus-fcpertminus)       &           ierr .EQ. 0 ) then
                      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  
                     
                      tmpplot1(itilepos,jtilepos,layer,itile,jtile) =  
      &                    gfd  
                      tmpplot2(itilepos,jtilepos,layer,itile,jtile) =  
      &                    ratio_ad  
                      tmpplot3(itilepos,jtilepos,layer,itile,jtile) =  
      &                    ratio_ftl  
362    
363                   if ( grdchk_eps .eq. 0. ) then
364                      gfd = (fcpertplus-fcpertminus)
365                   else
366                      gfd = (fcpertplus-fcpertminus)
367         &                 /(grdchk_epsfac*grdchk_eps)
368                   endif
369                                            
370                       fcrmem      ( ichknum ) = fcref                 if ( adxxmemo .eq. 0. ) then
371                       fcppmem     ( ichknum ) = fcpertplus                    ratio_ad = abs( adxxmemo - gfd )
                      fcpmmem     ( ichknum ) = fcpertminus  
                      xxmemref    ( ichknum ) = xxmemo_ref  
                      xxmempert   ( ichknum ) = xxmemo_pert  
                      gfdmem      ( ichknum ) = gfd  
                      adxxmem     ( ichknum ) = adxxmemo  
                      ftlxxmem    ( ichknum ) = ftlxxmemo  
                      ratioadmem  ( ichknum ) = ratio_ad  
                      ratioftlmem ( ichknum ) = ratio_ftl  
   
                      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 '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',  
      &                    ichknum,itilepos,jtilepos,layer,  
      &                    fcref, fcpertplus, fcpertminus  
                      print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',  
      &                    ichknum,ichkmem(ichknum),  
      &                    icompmem(ichknum),itestmem(ichknum),  
      &                    adxxmemo, gfd, ratio_ad  
                      print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',  
      &                    ichknum,ichkmem(ichknum),  
      &                    icompmem(ichknum),itestmem(ichknum),  
      &                    ftlxxmemo, gfd, ratio_ftl  
 cph)  
                   else  
 c  
                      print *, 'ph-grd 3 -------------------------------'  
                      print *, 'ph-grd 3 : ierr = ', ierr,  
      &                                 ', icomp = ', icomp  
                   endif  
372                 else                 else
373                    ierr_grdchk = -1                    ratio_ad = 1. - gfd/adxxmemo
374                 endif                 endif
               
             enddo  
          endif  
375    
376  c--   Everyone has to wait for the component to be reset.                 if ( ftlxxmemo .eq. 0. ) then
377           _BARRIER                    ratio_ftl = abs( ftlxxmemo - gfd )
378                   else
379                      ratio_ftl = 1. - gfd/ftlxxmemo
380                   endif
381                      
382                   tmpplot1(itilepos,jtilepos,layer,itile,jtile)
383         &              = gfd
384                   tmpplot2(itilepos,jtilepos,layer,itile,jtile)
385         &              = ratio_ad
386                   tmpplot3(itilepos,jtilepos,layer,itile,jtile)
387         &              = ratio_ftl
388    
389                   fcrmem      ( ichknum ) = fcref
390                   fcppmem     ( ichknum ) = fcpertplus
391                   fcpmmem     ( ichknum ) = fcpertminus
392                   xxmemref    ( ichknum ) = xxmemo_ref
393                   xxmempert   ( ichknum ) = xxmemo_pert
394                   gfdmem      ( ichknum ) = gfd
395                   adxxmem     ( ichknum ) = adxxmemo
396                   ftlxxmem    ( ichknum ) = ftlxxmemo
397                   ratioadmem  ( ichknum ) = ratio_ad
398                   ratioftlmem ( ichknum ) = ratio_ftl
399    
400                   irecmem   ( ichknum ) = icvrec
401                   bimem     ( ichknum ) = itile
402                   bjmem     ( ichknum ) = jtile
403                   ilocmem   ( ichknum ) = itilepos
404                   jlocmem   ( ichknum ) = jtilepos
405                   klocmem   ( ichknum ) = layer
406                   iobcsmem  ( ichknum ) = obcspos
407                   icompmem  ( ichknum ) = icomp
408                   ichkmem   ( ichknum ) = ichknum
409                   itestmem  ( ichknum ) = itest
410                   ierrmem   ( ichknum ) = ierr
411                      
412                   print *, 'grad-res -------------------------------'
413                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
414         &              myprocid,ichknum,itilepos,jtilepos,layer,
415         &              fcref, fcpertplus, fcpertminus
416    #ifdef ALLOW_TANGENTLINEAR_RUN
417                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
418         &              myprocid,ichknum,ichkmem(ichknum),
419         &              icompmem(ichknum),itestmem(ichknum),
420         &              ftlxxmemo, gfd, ratio_ftl
421                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
422         &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
423                   CALL PRINT_MESSAGE
424         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
425    #else
426                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
427         &              myprocid,ichknum,ichkmem(ichknum),
428         &              icompmem(ichknum),itestmem(ichknum),
429         &              adxxmemo, gfd, ratio_ad
430                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
431         &              'precision_grdchk_result ADM ', fcref, adxxmemo
432                   CALL PRINT_MESSAGE
433         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
434    #endif
435    
436                endif
437    
438                print *, 'ph-grd  ierr ---------------------------'
439                print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
440         &           ', ichknum = ', ichknum
441    
442                _BARRIER
443    
444    c-- else of if ( ichknum ...
445             else
446                ierr_grdchk = -1
447    
448    c-- end of if ( ichknum ...
449             endif
450    
451    c-- end of do icomp = ...
452        enddo        enddo
453    
454        CALL WRITE_REC_XYZ_RL( 'grd_findiff'   , tmpplot1, 1, 0, myThid)        if ( myProcId .EQ. grdchkwhichproc ) then
455        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)           CALL WRITE_REC_XYZ_RL(
456        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
457             CALL WRITE_REC_XYZ_RL(
458         &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
459             CALL WRITE_REC_XYZ_RL(
460         &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
461          endif
462    
463    c--   Everyone has to wait for the component to be reset.
464          _BARRIER
465    
466  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
467        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
468    
469  #endif /* ALLOW_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
470    
471        end        end
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22