/[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.5 by heimbach, Sat Jul 13 02:55:58 2002 UTC revision 1.11 by heimbach, Mon Oct 27 22:32:55 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 106  c     == local variables == Line 115  c     == local variables ==
115        _RL     gfd        _RL     gfd
116        _RL     fcref        _RL     fcref
117        _RL     fcpertplus, fcpertminus        _RL     fcpertplus, fcpertminus
118        _RL     ratio        _RL     ratio_ad
119          _RL     ratio_ftl
120        _RL     xxmemo_ref        _RL     xxmemo_ref
121        _RL     xxmemo_pert        _RL     xxmemo_pert
122        _RL     adxxmemo        _RL     adxxmemo
123          _RL     ftlxxmemo
124          _RL     localEps
125          _RL     grdchk_epsfac
126    
 cph(  
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  cph)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
   
 Cml(  
       _RL     grdchk_epsfac  
 Cml)  
130    
131  c     == end of interface ==  c     == end of interface ==
132  CEOP  CEOP
# Line 142  cph   Gradient via adjoint has already b Line 150  cph   Gradient via adjoint has already b
150  cph   and so has unperturbed cost function,  cph   and so has unperturbed cost function,
151  cph   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
152    
       fcref = fc  
153        ierr_grdchk = 0        ierr_grdchk = 0
154    cphadmtlm(
155          fcref = fc
156    cphadmtlm      fcref = objf_state_final(45,4,1,1)
157    cphadmtlm)
158    
159          print *, 'ph-check fcref = ', fcref
160    
 cph(  
161        do bj = jtlo, jthi        do bj = jtlo, jthi
162           do bi = itlo, ithi           do bi = itlo, ithi
163              do k = 1, nr              do k = 1, nr
# Line 153  cph( Line 165  cph(
165                    do i = imin, imax                    do i = imin, imax
166                       tmpplot1(i,j,k,bi,bj) = 0. _d 0                       tmpplot1(i,j,k,bi,bj) = 0. _d 0
167                       tmpplot2(i,j,k,bi,bj) = 0. _d 0                       tmpplot2(i,j,k,bi,bj) = 0. _d 0
168                         tmpplot3(i,j,k,bi,bj) = 0. _d 0
169                    end do                    end do
170                 end do                 end do
171              end do              end do
172           end do           end do
173        end do        end do
 cph)  
174    
175        if ( useCentralDiff ) then        if ( useCentralDiff ) then
176           grdchk_epsfac = 2. _d 0           grdchk_epsfac = 2. _d 0
# Line 166  cph) Line 178  cph)
178           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
179        end if        end if
180    
181          print *, 'ph-grd 3 -------------------------------'
182          print ('(2a)'),
183         &     ' ph-grd 3  proc    #    i    j    k            fc ref',
184         &     '        fc + eps        fc - eps'
185    #ifdef ALLOW_TANGENTLINEAR_RUN
186          print ('(2a)'),
187         &     ' ph-grd 3  proc    #    i    j    k          tlm grad',
188         &     '         fd grad      1 - fd/tlm'
189    #else
190          print ('(2a)'),
191         &     ' ph-grd 3  proc    #    i    j    k          adj grad',
192         &     '         fd grad      1 - fd/adj'
193    #endif
194    
195  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
196  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
197  c--   gradient checks.  c--   gradient checks.
       do jprocs = 1,numberOfProcs  
          proc_grdchk = jprocs - 1  
   
          if ( myProcId .eq. proc_grdchk ) then  
198    
199              do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
200    
201                 ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
202    
203                 if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
204    
205  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
206                    call grdchk_loc( icomp, ichknum,              if ( myProcId .EQ. grdchkwhichproc ) then
207       &                 icvrec, itile, jtile, layer,                 call grdchk_loc( icomp, ichknum,
208       &                 itilepos, jtilepos, itest, ierr,       &              icvrec, itile, jtile, layer, obcspos,
209       &                 mythid )       &              itilepos, jtilepos, itest, ierr,
210         &              mythid )
211                endif
212                _BARRIER
213                                
214                    if ( ierr .eq. 0 ) then  c******************************************************
215    c--   (A): get gradient component calculated via adjoint
216  c--   positive perturbation  c******************************************************
                      grdchk_eps = abs(grdchk_eps)  
217    
218  c--   get gradient component calculated via adjoint  c--   get gradient component calculated via adjoint
219                       call grdchk_getadxx( icvrec,              if ( myProcId .EQ. grdchkwhichproc .AND.
220       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
221       &                    itilepos, jtilepos,                 call grdchk_getadxx( icvrec,
222       &                    adxxmemo, mythid )       &              itile, jtile, layer,
223                       _BARRIER       &              itilepos, jtilepos,
224         &              adxxmemo, mythid )
225                endif
226                _BARRIER
227    
228    #ifdef ALLOW_TANGENTLINEAR_RUN
229    c******************************************************
230    c--   (B): Get gradient component g_fc from tangent linear run:
231    c******************************************************
232    c--
233    c--   1. perturb control vector component: xx(i)=1.
234    
235                if ( myProcId .EQ. grdchkwhichproc .AND.
236         &           ierr .EQ. 0 ) then
237                   localEps = 1. _d 0
238                   call grdchk_getxx( icvrec, TANGENT_SIMULATION,
239         &              itile, jtile, layer,
240         &              itilepos, jtilepos,
241         &              xxmemo_ref, xxmemo_pert, localEps,
242         &              mythid )
243                endif
244                _BARRIER
245    
246    c--
247    c--   2. perform tangent linear run
248                mytime = starttime
249                myiter = niter0
250    cphadmtlm(
251                g_fc = 0.
252    cphadmtlm            do j=1,sny
253    cphadmtlm               do i=1,snx
254    cphadmtlm                  g_objf_state_final(i,j,1,1) = 0.
255    cphadmtlm               enddo
256    cphadmtlm            enddo
257    cphadmtlm)
258                call g_the_main_loop( mytime, myiter, mythid )
259    cphadmtlm(
260                ftlxxmemo = g_fc
261    cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1)
262    cphadmtlm)
263                _BARRIER
264    c--
265    c--   3. reset control vector
266                if ( myProcId .EQ. grdchkwhichproc .AND.
267         &           ierr .EQ. 0 ) then
268                   call grdchk_setxx( icvrec, TANGENT_SIMULATION,
269         &              itile, jtile, layer,
270         &              itilepos, jtilepos,
271         &              xxmemo_ref, mythid )
272                endif
273                _BARRIER
274    
275    #endif /* ALLOW_TANGENTLINEAR_RUN */
276    
277    
278    c******************************************************
279    c--   (C): Get gradient via finite difference perturbation
280    c******************************************************
281    
282  c--   get control vector component from file  c--   get control vector component from file
283  c--   perturb it and write back to file  c--   perturb it and write back to file
284                       call grdchk_getxx( icvrec,  c--   positive perturbation
285       &                    itile, jtile, layer,              localEps = abs(grdchk_eps)
286       &                    itilepos, jtilepos,              if ( myProcId .EQ. grdchkwhichproc .AND.
287       &                    xxmemo_ref, xxmemo_pert, mythid )       &           ierr .EQ. 0 ) then
288                       _BARRIER                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
289         &              itile, jtile, layer,
290         &              itilepos, jtilepos,
291         &              xxmemo_ref, xxmemo_pert, localEps,
292         &              mythid )
293                endif
294                _BARRIER
295                                            
296  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
297                       mytime = starttime              mytime = starttime
298                       myiter = niter0              myiter = niter0
299                       call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
300                       fcpertplus = fc  cphadmtlm(
301                fcpertplus = fc
302    cphadmtlm            fcpertplus = objf_state_final(45,4,1,1)
303    cphadmtlm)
304                _BARRIER
305                                        
306  c--   Reset control vector.  c--   Reset control vector.
307                       call grdchk_setxx( icvrec,              if ( myProcId .EQ. grdchkwhichproc .AND.
308       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
309       &                    itilepos, jtilepos,                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
310       &                    xxmemo_ref, mythid )       &              itile, jtile, layer,
311                       _BARRIER       &              itilepos, jtilepos,
312         &              xxmemo_ref, mythid )
313                endif
314                _BARRIER
315    
316                       fcpertminus = fcref              fcpertminus = fcref
317    
318                       if ( useCentralDiff ) then              if ( useCentralDiff ) then
   
 c--   repeat the proceedure for a negative perturbation  
                         grdchk_eps = - abs(grdchk_eps)  
319    
320  c--   get control vector component from file  c--   get control vector component from file
321  c--   perturb it and write back to file  c--   perturb it and write back to file
322                          call grdchk_getxx( icvrec,  c--   repeat the proceedure for a negative perturbation
323       &                    itile, jtile, layer,                 if ( myProcId .EQ. grdchkwhichproc .AND.
324       &                    itilepos, jtilepos,       &           ierr .EQ. 0 ) then
325       &                    xxmemo_ref, xxmemo_pert, mythid )                    localEps = - abs(grdchk_eps)
326                          _BARRIER                    call grdchk_getxx( icvrec, FORWARD_SIMULATION,
327         &                 itile, jtile, layer,
328         &                 itilepos, jtilepos,
329         &                 xxmemo_ref, xxmemo_pert, localEps,
330         &                 mythid )
331                   endif
332                   _BARRIER
333                                            
334  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
335                          mytime = starttime                 mytime = starttime
336                          myiter = niter0                 myiter = niter0
337                          call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
338                          fcpertminus = fc                 _BARRIER
339                   fcpertminus = fc
340                                        
341  c--   Reset control vector.  c--   Reset control vector.
342                          call grdchk_setxx( icvrec,                 if ( myProcId .EQ. grdchkwhichproc .AND.
343       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
344       &                    itilepos, jtilepos,                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,
345       &                    xxmemo_ref, mythid )       &                 itile, jtile, layer,
346                          _BARRIER       &                 itilepos, jtilepos,
347         &                 xxmemo_ref, mythid )
348                   endif
349                   _BARRIER
350    
351  c     reset grdchk_eps to a positive value  c-- end of if useCentralDiff ...
352                          grdchk_eps = abs(grdchk_eps)              end if
353    
354                       end if  c******************************************************
355  c--  c--   (D): calculate relative differences between gradients
356                       if ( grdchk_eps .eq. 0. ) then  c******************************************************
                         gfd = (fcpertplus-fcpertminus)  
                      else  
                         gfd = (fcpertplus-fcpertminus)  
      &                       /(grdchk_epsfac*grdchk_eps)  
                      endif  
                       
                      if ( adxxmemo .eq. 0. ) then  
                         ratio = abs( adxxmemo - gfd )  
                      else  
                         ratio = 1. - gfd/adxxmemo  
                      endif  
                     
 cph(  
                      tmpplot1(itilepos,jtilepos,layer,itile,jtile) =  
      &                    gfd  
                      tmpplot2(itilepos,jtilepos,layer,itile,jtile) =  
      &                    ratio  
 cph)  
357    
358                if ( myProcId .EQ. grdchkwhichproc .AND.
359         &           ierr .EQ. 0 ) then
360    
361                   if ( grdchk_eps .eq. 0. ) then
362                      gfd = (fcpertplus-fcpertminus)
363                   else
364                      gfd = (fcpertplus-fcpertminus)
365         &                 /(grdchk_epsfac*grdchk_eps)
366                   endif
367                                            
368                       fcrmem    ( ichknum ) = fcref                 if ( adxxmemo .eq. 0. ) then
369                       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  
                      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 '(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  
 cph)  
                   else  
 c  
                      print *, 'ph-grd 3 -------------------------------'  
                      print *, 'ph-grd 3 : ierr = ', ierr,  
      &                                 ', icomp = ', icomp  
                   endif  
370                 else                 else
371                    ierr_grdchk = -1                    ratio_ad = 1. - gfd/adxxmemo
372                 endif                 endif
               
             enddo  
          endif  
373    
374  c--   Everyone has to wait for the component to be reset.                 if ( ftlxxmemo .eq. 0. ) then
375           _BARRIER                    ratio_ftl = abs( ftlxxmemo - gfd )
376                   else
377                      ratio_ftl = 1. - gfd/ftlxxmemo
378                   endif
379                      
380                   tmpplot1(itilepos,jtilepos,layer,itile,jtile)
381         &              = gfd
382                   tmpplot2(itilepos,jtilepos,layer,itile,jtile)
383         &              = ratio_ad
384                   tmpplot3(itilepos,jtilepos,layer,itile,jtile)
385         &              = ratio_ftl
386    
387                   fcrmem      ( ichknum ) = fcref
388                   fcppmem     ( ichknum ) = fcpertplus
389                   fcpmmem     ( ichknum ) = fcpertminus
390                   xxmemref    ( ichknum ) = xxmemo_ref
391                   xxmempert   ( ichknum ) = xxmemo_pert
392                   gfdmem      ( ichknum ) = gfd
393                   adxxmem     ( ichknum ) = adxxmemo
394                   ftlxxmem    ( ichknum ) = ftlxxmemo
395                   ratioadmem  ( ichknum ) = ratio_ad
396                   ratioftlmem ( ichknum ) = ratio_ftl
397    
398                   irecmem   ( ichknum ) = icvrec
399                   bimem     ( ichknum ) = itile
400                   bjmem     ( ichknum ) = jtile
401                   ilocmem   ( ichknum ) = itilepos
402                   jlocmem   ( ichknum ) = jtilepos
403                   klocmem   ( ichknum ) = layer
404                   iobcsmem  ( ichknum ) = obcspos
405                   icompmem  ( ichknum ) = icomp
406                   ichkmem   ( ichknum ) = ichknum
407                   itestmem  ( ichknum ) = itest
408                   ierrmem   ( ichknum ) = ierr
409                      
410                   print *, 'ph-grd 3 -------------------------------'
411                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
412         &              myprocid,ichknum,itilepos,jtilepos,layer,
413         &              fcref, fcpertplus, fcpertminus
414    #ifdef ALLOW_TANGENTLINEAR_RUN
415                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
416         &              myprocid,ichknum,ichkmem(ichknum),
417         &              icompmem(ichknum),itestmem(ichknum),
418         &              ftlxxmemo, gfd, ratio_ftl
419    #else
420                   print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
421         &              myprocid,ichknum,ichkmem(ichknum),
422         &              icompmem(ichknum),itestmem(ichknum),
423         &              adxxmemo, gfd, ratio_ad
424    #endif
425    
426                endif
427    
428                print *, 'ph-grd  ierr ---------------------------'
429                print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
430         &           ', ichknum = ', ichknum
431    
432                _BARRIER
433    
434    c-- else of if ( ichknum ...
435             else
436                ierr_grdchk = -1
437    
438    c-- end of if ( ichknum ...
439             endif
440    
441    c-- end of do icomp = ...
442        enddo        enddo
443    
444  cph(        if ( myProcId .EQ. grdchkwhichproc ) then
445        CALL WRITE_REC_XYZ_RL( 'grd_findiff' , tmpplot1, 1, 0, myThid)           CALL WRITE_REC_XYZ_RL(
446        CALL WRITE_REC_XYZ_RL( 'grd_ratio' , tmpplot2, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
447  cph)           CALL WRITE_REC_XYZ_RL(
448         &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
449             CALL WRITE_REC_XYZ_RL(
450         &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
451          endif
452    
453    c--   Everyone has to wait for the component to be reset.
454          _BARRIER
455    
456  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
457        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
458    
459  #endif /* ALLOW_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
460    
461        end        end
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22