/[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.15 by heimbach, Tue Oct 26 20:10:25 2004 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,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  
200    
201           if ( myProcId .eq. proc_grdchk ) then        if ( nbeg .EQ. 0 ) call grdchk_get_position( mythid )
202    
203              do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
204    
205                 ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
206    
207                 if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
208    
209  c--   Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
210                    call grdchk_loc( icomp, ichknum,              if ( myProcId .EQ. grdchkwhichproc ) then
211       &                 icvrec, itile, jtile, layer,                 call grdchk_loc( icomp, ichknum,
212       &                 itilepos, jtilepos, itest, ierr,       &              icvrec, itile, jtile, layer, obcspos,
213       &                 mythid )       &              itilepos, jtilepos, itest, ierr,
214         &              mythid )
215                endif
216                _BARRIER
217                                
                   if ( ierr .eq. 0 ) then  
   
218  c******************************************************  c******************************************************
219  c--   (A): get gradient component calculated via adjoint  c--   (A): get gradient component calculated via adjoint
220  c******************************************************  c******************************************************
221                       call grdchk_getadxx( icvrec,  
222       &                    itile, jtile, layer,  c--   get gradient component calculated via adjoint
223       &                    itilepos, jtilepos,              if ( myProcId .EQ. grdchkwhichproc .AND.
224       &                    adxxmemo, mythid )       &           ierr .EQ. 0 ) then
225                       _BARRIER                 call grdchk_getadxx( icvrec,
226         &              itile, jtile, layer,
227         &              itilepos, jtilepos,
228         &              adxxmemo, mythid )
229                endif
230                _BARRIER
231    
232  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
233  c******************************************************  c******************************************************
# Line 207  c*************************************** Line 236  c***************************************
236  c--  c--
237  c--   1. perturb control vector component: xx(i)=1.  c--   1. perturb control vector component: xx(i)=1.
238    
239                       localEps = 1.              if ( myProcId .EQ. grdchkwhichproc .AND.
240                       call grdchk_getxx( icvrec, TANGENT_SIMULATION,       &           ierr .EQ. 0 ) then
241       &                    itile, jtile, layer,                 localEps = 1. _d 0
242       &                    itilepos, jtilepos,                 call grdchk_getxx( icvrec, TANGENT_SIMULATION,
243       &                    xxmemo_ref, xxmemo_pert, localEps,       &              itile, jtile, layer,
244       &                    mythid )       &              itilepos, jtilepos,
245                       _BARRIER       &              xxmemo_ref, xxmemo_pert, localEps,
246         &              mythid )
247                endif
248                _BARRIER
249    
250  c--  c--
251  c--   2. perform tangent linear run  c--   2. perform tangent linear run
252                       mytime = starttime              mytime = starttime
253                       myiter = niter0              myiter = niter0
254                       g_fc = 0.  cphadmtlm(
255                       call g_the_main_loop( mytime, myiter, mythid )              g_fc = 0.
256                       ftlxxmemo = g_fc  cphadmtlm            do j=1,sny
257    cphadmtlm               do i=1,snx
258    cphadmtlm                  g_objf_state_final(i,j,1,1,1) = 0.
259    cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 0.
260    cphadmtlm               enddo
261    cphadmtlm            enddo
262    cphadmtlm)
263                call g_the_main_loop( mytime, myiter, mythid )
264    cphadmtlm(
265                ftlxxmemo = g_fc
266    cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1,1)
267    cphadmtlm)
268                _BARRIER
269  c--  c--
270  c--   3. reset control vector  c--   3. reset control vector
271                       call grdchk_setxx( icvrec, TANGENT_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
272       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
273       &                    itilepos, jtilepos,                 call grdchk_setxx( icvrec, TANGENT_SIMULATION,
274       &                    xxmemo_ref, mythid )       &              itile, jtile, layer,
275                       _BARRIER       &              itilepos, jtilepos,
276  #endif       &              xxmemo_ref, mythid )
277                endif
278                _BARRIER
279    
280    #endif /* ALLOW_TANGENTLINEAR_RUN */
281    
282    
283  c******************************************************  c******************************************************
284  c--   (C): Get gradient via finite difference perturbation  c--   (C): Get gradient via finite difference perturbation
285  c******************************************************  c******************************************************
286    
287  c--   get control vector component from file  c--   get control vector component from file
288  c--   perturb it and write back to file:  c--   perturb it and write back to file
289  c--   positive perturbation  c--   positive perturbation
290                       localEps = abs(grdchk_eps)              localEps = abs(grdchk_eps)
291                       call grdchk_getxx( icvrec, FORWARD_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
292       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
293       &                    itilepos, jtilepos,                 call grdchk_getxx( icvrec, FORWARD_SIMULATION,
294       &                    xxmemo_ref, xxmemo_pert, localEps,       &              itile, jtile, layer,
295       &                    mythid )       &              itilepos, jtilepos,
296                       _BARRIER       &              xxmemo_ref, xxmemo_pert, localEps,
297         &              mythid )
298                endif
299                _BARRIER
300                                            
301  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
302                       mytime = starttime              mytime = starttime
303                       myiter = niter0              myiter = niter0
304                       call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
305                       fcpertplus = fc  cphadmtlm(
306                fcpertplus = fc
307    cphadmtlm            fcpertplus = objf_state_final(45,4,1,1,1)
308    cphadmtlm)
309                print *, 'ph-check fcpertplus = ', fcpertplus
310                _BARRIER
311                                        
312  c--   Reset control vector.  c--   Reset control vector.
313                       call grdchk_setxx( icvrec, FORWARD_SIMULATION,              if ( myProcId .EQ. grdchkwhichproc .AND.
314       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
315       &                    itilepos, jtilepos,                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,
316       &                    xxmemo_ref, mythid )       &              itile, jtile, layer,
317                       _BARRIER       &              itilepos, jtilepos,
318         &              xxmemo_ref, mythid )
319                endif
320                _BARRIER
321    
322                       fcpertminus = fcref              fcpertminus = fcref
323                print *, 'ph-check fcpertminus = ', fcpertminus
324    
325                       if ( useCentralDiff ) then              if ( useCentralDiff ) then
326    
327  c--   get control vector component from file  c--   get control vector component from file
328  c--   perturb it and write back to file:  c--   perturb it and write back to file
329  c--   repeat the proceedure for a negative perturbation  c--   repeat the proceedure for a negative perturbation
330                          localEps = - abs(grdchk_eps)                 if ( myProcId .EQ. grdchkwhichproc .AND.
331                          call grdchk_getxx( icvrec, FORWARD_SIMULATION,       &           ierr .EQ. 0 ) then
332       &                    itile, jtile, layer,                    localEps = - abs(grdchk_eps)
333       &                    itilepos, jtilepos,                    call grdchk_getxx( icvrec, FORWARD_SIMULATION,
334       &                    xxmemo_ref, xxmemo_pert, localEps,       &                 itile, jtile, layer,
335       &                    mythid )       &                 itilepos, jtilepos,
336                          _BARRIER       &                 xxmemo_ref, xxmemo_pert, localEps,
337         &                 mythid )
338                   endif
339                   _BARRIER
340                                            
341  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
342                          mytime = starttime                 mytime = starttime
343                          myiter = niter0                 myiter = niter0
344                          call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
345                          fcpertminus = fc                 _BARRIER
346                   fcpertminus = fc
347                                        
348  c--   Reset control vector.  c--   Reset control vector.
349                          call grdchk_setxx( icvrec, FORWARD_SIMULATION,                 if ( myProcId .EQ. grdchkwhichproc .AND.
350       &                    itile, jtile, layer,       &           ierr .EQ. 0 ) then
351       &                    itilepos, jtilepos,                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,
352       &                    xxmemo_ref, mythid )       &                 itile, jtile, layer,
353                          _BARRIER       &                 itilepos, jtilepos,
354         &                 xxmemo_ref, mythid )
355                   endif
356                   _BARRIER
357    
358    c-- end of if useCentralDiff ...
359                end if
360    
                      end if  
 c--  
361  c******************************************************  c******************************************************
362  c--   (D): calculate relative differences between gradients  c--   (D): calculate relative differences between gradients
363  c******************************************************  c******************************************************
364    
365                       if ( grdchk_eps .eq. 0. ) then              if ( myProcId .EQ. grdchkwhichproc .AND.
366                          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  
367    
368                   if ( grdchk_eps .eq. 0. ) then
369                      gfd = (fcpertplus-fcpertminus)
370                   else
371                      gfd = (fcpertplus-fcpertminus)
372         &                 /(grdchk_epsfac*grdchk_eps)
373                   endif
374                                            
375                       fcrmem      ( ichknum ) = fcref                 if ( adxxmemo .eq. 0. ) then
376                       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  
377                 else                 else
378                    ierr_grdchk = -1                    ratio_ad = 1. - gfd/adxxmemo
379                 endif                 endif
               
             enddo  
          endif  
380    
381  c--   Everyone has to wait for the component to be reset.                 if ( ftlxxmemo .eq. 0. ) then
382           _BARRIER                    ratio_ftl = abs( ftlxxmemo - gfd )
383                   else
384                      ratio_ftl = 1. - gfd/ftlxxmemo
385                   endif
386                      
387                   tmpplot1(itilepos,jtilepos,layer,itile,jtile)
388         &              = gfd
389                   tmpplot2(itilepos,jtilepos,layer,itile,jtile)
390         &              = ratio_ad
391                   tmpplot3(itilepos,jtilepos,layer,itile,jtile)
392         &              = ratio_ftl
393    
394                   fcrmem      ( ichknum ) = fcref
395                   fcppmem     ( ichknum ) = fcpertplus
396                   fcpmmem     ( ichknum ) = fcpertminus
397                   xxmemref    ( ichknum ) = xxmemo_ref
398                   xxmempert   ( ichknum ) = xxmemo_pert
399                   gfdmem      ( ichknum ) = gfd
400                   adxxmem     ( ichknum ) = adxxmemo
401                   ftlxxmem    ( ichknum ) = ftlxxmemo
402                   ratioadmem  ( ichknum ) = ratio_ad
403                   ratioftlmem ( ichknum ) = ratio_ftl
404    
405                   irecmem   ( ichknum ) = icvrec
406                   bimem     ( ichknum ) = itile
407                   bjmem     ( ichknum ) = jtile
408                   ilocmem   ( ichknum ) = itilepos
409                   jlocmem   ( ichknum ) = jtilepos
410                   klocmem   ( ichknum ) = layer
411                   iobcsmem  ( ichknum ) = obcspos
412                   icompmem  ( ichknum ) = icomp
413                   ichkmem   ( ichknum ) = ichknum
414                   itestmem  ( ichknum ) = itest
415                   ierrmem   ( ichknum ) = ierr
416                      
417                   print *, 'grad-res -------------------------------'
418                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
419         &              myprocid,ichknum,itilepos,jtilepos,layer,
420         &              fcref, fcpertplus, fcpertminus
421    #ifdef ALLOW_TANGENTLINEAR_RUN
422                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
423         &              myprocid,ichknum,ichkmem(ichknum),
424         &              icompmem(ichknum),itestmem(ichknum),
425         &              ftlxxmemo, gfd, ratio_ftl
426                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
427         &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
428                   CALL PRINT_MESSAGE
429         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
430    #else
431                   print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',
432         &              myprocid,ichknum,ichkmem(ichknum),
433         &              icompmem(ichknum),itestmem(ichknum),
434         &              adxxmemo, gfd, ratio_ad
435                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
436         &              'precision_grdchk_result ADM ', fcref, adxxmemo
437                   CALL PRINT_MESSAGE
438         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
439    #endif
440    
441                endif
442    
443                print *, 'ph-grd  ierr ---------------------------'
444                print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
445         &           ', ichknum = ', ichknum
446    
447                _BARRIER
448    
449    c-- else of if ( ichknum ...
450             else
451                ierr_grdchk = -1
452    
453    c-- end of if ( ichknum ...
454             endif
455    
456    c-- end of do icomp = ...
457        enddo        enddo
458    
459        CALL WRITE_REC_XYZ_RL( 'grd_findiff'   , tmpplot1, 1, 0, myThid)        if ( myProcId .EQ. grdchkwhichproc ) then
460        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)           CALL WRITE_REC_XYZ_RL(
461        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)       &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
462             CALL WRITE_REC_XYZ_RL(
463         &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
464             CALL WRITE_REC_XYZ_RL(
465         &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
466          endif
467    
468    c--   Everyone has to wait for the component to be reset.
469          _BARRIER
470    
471  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
472        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
473    
474  #endif /* ALLOW_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
475    
476        end        end
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22