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